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PREFACE 


The  investigation  reported  herein  was  conducted  for  and  sponsored 
by  the  U.S.  Army  Engineer  Waterways  Experiment  Station  (WES,  Contract 
No.  DACW39-80-C-0129).  The  numerical  aspects  of  the  investigation 
served  as  the  basis  for  the  Ph.O.  thesis  of  Mr.  John  Vadnal.  The 
companion,  experimental  investigation  has  been  reported  by  Odgaard  and 
Kennedy 3.  The  numerical  model  developed  in  this  phase  of  the 

investigation  analyzes  and  predicts  flow  and  sediment-transport 
distributions  in  alluvial-channel  bends.  The  authors  wish  to 
acknowledge  their  gratitude  to  Mr.  Steve  Maynord  of  WES  for  his 
continuing  encouragement  and  assistance  during  the  course  of  this  study. 
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CONVERSION  FACTORS,  U.S.  CUSTOMARY  TO  METRIC(SI) 

UNITS  OF  MEASUREMENT 

U.S.  customary  units  of  measurement  used  in  this  report  can  hP 
converted  to  metric  (SI)  units  as  follows:  P  e 


Multiply 


To  Obtain 


inches 

feet 

feet  per  second 

cubic  feet  per  second 

miles  (U.S.  statute) 

pounds  (mass) 

pounds  (force)  per  square 

Inch 

tons  (short)  per  foot 
per  day 


0.0254 

0.3048 

0.3048 

0.02831685 

1.6093 

0.4535924 

6894.757 

2.9763 


meters 

meters 

meters  per  second 
cubic  meters  per  second 
ki 1 ometers 
kilograms 

pascals 

tons  per  meter  per  day 


A  NUMERICAL  MODEL  FOR  FLOW  AND  SEDIMENT  TRANSPORT  IN 
ALLUVIAL-RIVER  BENDS 


PART  I:  INTRODUCTION 
Background 

1.  Two  of  the  most  striking  and  intriguing  features  of  natural 
alluvial  streams  are  their  tendency  to  meander,  and  the  downstream 
migration  of  the  meanders.  In  addition  to  being  a  fascinating  natural 
phenomenon  and  posing  some  of  the  most  nettling  problems  in  the  whole  of 
river  mechanics,  river  meandering,  and  in  particular  the  bank  erosion 
attendant  to  the  growth  and  migration  of  the  meander  loops,  has  become  a 
major  International  problem.  According  to  the  final  report  on  work 
conducted  under  the  Streambank  Erosion  Control  Evaluation  and 
Demonstration  Act  of  1974*  (Section  32,  Public  Law  32-251,  submitted  In 
December  1981),  approximately  142,000  bank-miles  of  streams  and 
waterways  are  in  need  of  erosion  protection.  The  cost  to  arrest  or 
control  this  erosion  by  means  of  conventional  bank -protect! on  methods 
currently  available  is  estimated  to  be  In  excess  of  $1  billion 
annually.  For  the  Upper-Mi ssissippi  River  basin  alone,  the  cost  was 
estimated  to  be  in  excess  of  $21  million  annually.  These  figures  exceed 
the  benefits  derived  by  a  large  margin,  thereby  rendering  many  of  the 
bank-eroslon-control  projects  uneconomical  on  a  cost/benefit  basis.  As 
a  result,  most  bank-erosion  losses  continue  unabated.  Attempts  to  halt 
the  erosion  are  often  limited  to  piecemeal  protection  along  Isolated 
bank  reaches,  at  public  or  private  facilities  on  streambanks,  or  at 
highway  crossings.  However,  as  such  facilities  increase  in  value  and  as 
the  consequences  of  failure  become  greater,  the  threshold  level  of 
acceptable  risk  becomes  smaller.  At  the  same  time,  traditional  channel  - 
stablllzatlon  measures  have  become  extremely  expensive  and  are  not 
acceptable  to  environmentalists  in  many  instances. 

2.  Nowhere  has  the  problem  come  to  sharper  focus  than  along  the 
Sacramento  River,  California.  The  Sacramento  River  Valley  contains  the 


Nation's  finest  (and  most  rapidly  disappearing)  agricultural  land. 
According  to  Brice^,  river-bank  erosion  along  the  unprotected  stretches 
of  the  approximately  200-mile-long  reach  of  the  Middle  and  Lower 
Sacramento  River  is  producing  an  average  annual  loss  of  nearly  two  acres 
of  farmland  per  mile.  Even  when  evaluated  against  current  inflated  land 
values*  traditional  means  of  bank  protection  (for  example  rock  riprap) 
are  so  expensive  that  they  cannot  be  justified  economically.  The 
problem  is  compounded  by  some  of  the  material  eroded  from  the  banks 
being  transported  to  the  dredged  navigation  channels  of  the  Lower 
Sacramento  River  system  and  San  Francisco  Bay.  Bank  protection  along 
the  upper  reaches  by  traditional  means  can  be  justified  economically 
only  if  it  can  be  demonstrated  that  the  reduced  erosion  will  result  in 
less  dredging  for  navi gati on -channel  maintenance.  Thus,  the  problem 
poses  two  general  questions:  (1)  Will  it  be  possible  to  develop 
alternative  bank-protection  measures  that  are  effective,  environmentally 
acceptable,  and  economically  justifiable  when  evaluated  against  land 
values  alone?;  and  (2)  Will  reduced  bank  erosion  upstream  be  reflected 
in  reduced  downstream  dredging  (and  how  much,  and  when),  or  is  the 
material  eroded  from  the  banks  being  deposited  at  other  locations  (e.g., 
on  point  bars)  along  the  river? 

3.  It  is  against  this  backdrop  that  the  Institute  of  Hydraulic 
Research  at  The  University  of  Iowa  entered  into  a  contract  with  the  Army 
Corps  of  Engineers,  Sacramento  District,  in  1980  to  conduct  an 
investigation  directed  at  developing  improved,  "unconventional"  bank- 
protection  methods  for  application  along  the  Sacramento  River.  It  was 
realized  that  the  investigation  should  also  include  laboratory  testing 
of  the  techniques  or  the  devices  proposed,  and  development  of  an 
analytical  model,  likely  a  computer-based  one,  for  routing  of  flow  and 
sediment  through  channel  bends.  Funds  for  conduct  of  the  laboratory 
investigation  and  development  of  the  routing  model  were  not  available 
from  the  Sacramento  District,  but  were  provided  by  the  Waterways 
Experiment  Station. 
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4.  A  report  describing  the  new  bank-protection  method  developed 
under  the  contract  with  the  Sacramento  District  (Contract  No.  DACW05-80- 
C-0083)  and  the  laboratory  testing  supported  jointly  by  the  Sacramento 
District  and  the  Waterways  Experiment  Station  (Contract  No.  DACW39-80-C- 
0129)  was  submitted  in  May  19823.  The  present  report  is  concerned 
solely  with  the  second  phase  of  the  WES-sponsored  project:  development 
of  a  numerical  model  for  analysis  and  prediction  of  flow  and  sediment- 
transport  characteristics  in  the  bends  of  meandering  alluvial  channels. 

Analytical  Strategy 

5.  The  point  of  departure  for  development  of  the  numerical  model 
is  the  analytical  work  reported  by  Falcon  and  Kennedy^.  The  manuscript 
of  this  paper  is  included  herewith  as  Appendix  A,  and  is  to  be 
considered  an  integral  part  of  this  report.  An  understanding  of  the 
analysis  presented  in  PART  II  requires  considerable  familiarity  with  the 
Falcon -Kennedy  analysis  described  in  Appendix  A. 

6.  The  principal  stumbling  blocks  encountered  in  the  analysis  of 
river-bend  flow  are  the  interdependency  of  the  bed  topography,  flow 
distribution,  and  sediment-discharge  distribution.  The  interaction  of 
the  vertically  nonuniform  distribution  of  streamwise  velocity  and 
channel  curvature  produces  a  secondary  flow  which  spirals  about  the 
channel -section  axis  and  moves  the  higher  velocity,  near-surface  fluid 
toward  the  outside  of  the  bend  and  the  near-bed  fluid  toward  the 
inside.  The  radially  inward  bed  shear-stress  transports  bed  sediment 
radially  inward  until  the  bed  becomes  inclined  such  that  the  radial- 
plane  shear  stress  is  balanced  by  the  component  of  the  moving  bed 
layer's  weight  radially  outward  along  the  bed.  The  resulting  warping  of 
the  channel  bed,  which  produces  larger  depths  near  the  outer  bank,  as 
shown  in  figure  1,  leads  to  a  redistr' bution  of  the  streamwise  flow,  and 
produces  much  larger  velocities,  boundary  shear  stress,  and  unit 
discharge  near  the  outer  bank.  This,  in  turn,  affects  the  lateral 
distribution  of  unit  sediment  discharge.  It  is  emphasized  that  the 
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secondary  flow  itself  has  a  relatively  minor  impact  on  the  distribution 
of  flow  and  sediment  transport  in  a  channel  bend.  It  is  the  bed  warping 
produced  by  the  secondary  flow  that  is  primarily  responsible  for  the 
redistributions  outlined  above. 

7.  In  the  case  of  fully  developed  flow  in  a  uniform  curved 
channel  ( i . e. »  one  of  whose  channel  axis  is  circular  in  plan  view),  the 
torque  generated  about  the  channel  axis  by  the  interaction  of  the 
velocity  profile  and  channel  curvature  is  balanced  almost  exclusively  by 
boundary  shear  stress.  In  the  case  of  channels  with  nonuniform  plan- 
form  curvature,  as  occurs  in  meandering  streams,  it  is  necessary  in  the 
calculation  of  the  secondary-f 1 ow  strength  to  include  the  nonuni fonni ty 
of  the  flux  of  moment -of -moment  urn  (or,  more  simply,  the  torsional 
inertia  of  the  flowing  fluid)  in  the  torque  balance.  The  inclusion  of 
the  inertial  effects  in  this  case  leads  to  a  phase  shift  between  local 
channel  curvature  and  local  secondary -current  strength  and  transverse 
bed  slope. 

8.  A  hallmark  of  the  numerical  flow  and  sediment  routing  model 
developed  in  PART  II  is  a  partial  uncoupling  of  the  secondary  flow  from 
the  distributions  of  primary  flow  and  sediment  transport.  The  principal 
steps  in  the  development  of  the  analytical  framework  for  the  numerical 
model  are  summarized  with  annotation  as  follows: 

i.  The  strength  of  the  secondary  flow  is  computed  for  any 

section  along  a  channel  of  nonuniform  curvature.  The 
integral -form  analysis  of  conservation  of  moment -of -momentum 
(or  torsion)  is  performed  for  a  channel  of  rectangular  cross 
section  with  depth  equal  to  that  of  the  same  flow  in  a 
straight  rectangular  channel  of  equal  slope. 

ii.  It  is  assumed  that  the  bed  topography  can  be  adequately 

represented  by  an  inclined  straight  line  passing  through  the 
mid-width  point  of  the  equivalent  rectangular  channel 


7 


utilized  in  step  i,  above.  It  is  further  argued  that  the 
transverse  slope  of  the  channel  bed  varies  linearly  with  the 
strength  of  the  secondary  current.  The  force  equilibrium  of 
the  bed  layer  Is  analyzed  to  relate  the  mean  bed  slope  to 
the  secondary-flow  strength. 

iii.  The  depth-integrated  differential  equations  of  continuity 
and  of  conservation  of  streamwise  momentum  then  are  employed 
to  calculate  the  velocity  field.  It  is  assumed  that  the 
secondary-flow  velocity  has  linear  vertical  distribution  at 
any  point  across  the  channel,  and  the  magnitu''  of  the 
secondary -flow  velocity  is  obtained  from  the  analysis 
presented  in  Appendix  A,  A  radial,  mass-shift  \  ocity  is 
also  included  to  account  for  the  movement  of  f  across 
the  channel  as  the  channel  curvature  and  trans  .  se  bed 
slope  change  and  even  reverse.  The  mass-shift  velocity  is 
assumed  to  be  uniformly  distributed  over  the  depth. 

iv.  The  lateral  distribution  of  unit  sediment  discharge  across 
the  channel  at  any  section  is  computed  on  the  basis  of  a 
power  law  using  the  local  flow  properties  computed  in  step 
iii. 

9.  The  analysis  is  limited  to  steady  flow  in  channels  with 
constant  width  and  centerline  streamwise  slope.  Utilization  of  the 
model  requires  an  estimate  of  the  stream's  friction  factor  using  some 
other  method,  such  as  that  of  Alam  and  Lovera^.  The  model  does  not 
allow  for  transverse  variations  of  local  friction  factor,  due  to  the 
lateral  variations  of  local  depth  and  velocity,  nor  does  it  compute 
sediment  discharge  by  size  fraction.  However,  the  computer  program  is 
structured  to  accommodate  these  features. 
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PART  II:  ANALYTICAL  MODEL 


Secondary  Flow  in  Rectangular  Channels  with  Nonuniform  Curvature 


10.  As  indicated  above,  the  secondary-flow  calculation  will  be 
made  for  a  rectangular  channel  that  is  "equivalent"  to  the  warped 
sections.  Because  the  analysis  is  of  the  integral  form,  and  considers 
only  the  streamwise  and  lateral  fluxes,  over  the  whole  channel  section, 
of  the  quantities  of  interest,  neglect  of  the  lateral  variations  of  the 
primary  and  secondary  currents  that  occur  in  a  warped  channel  is  not 
judged  to  be  a  major  limitation.  Support  for  this  conclusion  is 
provided  by  the  generally  good  agreement  between  measured  transverse  bed 
slopes  and  those  calculated  on  the  basis  of  the  radial  bed  shear  stress 
in  the  "equivalent"  rectangular  channel  (see  figure  4A). 

11.  The  control  volume  to  be  utilized  in  the  moment -of -moment urn 
analysis  is  depicted  in  figure  2,  and  the  coordinate  system  is  defined 
in  figure  3.  The  control  volume  can  be  envisioned  as  the  central  region 
(see  figure  1)  of  the  flow  as  it  existed  before  the  bed  became  warped  by 
the  secondary  current.  The  primary-flow  velocity  profile  will  be 
described  by  the  power  law, 


v.  _  n+1  1/n 

V  '  n  'd' 


(1) 


where,  in  addition  to  the  quantities  defined*  in  figures  1  and  2,  V  = 
depth -averaged  flow  velocity  and  n  =  reciprocal  of  the  power-law 
exponent.  The  secondary -flow  velocity  distribution  will  be  approximated 
by  a  linear  profile. 


(2) 


*For  convenience,  symbols  and  unusual  abbreviations  are  listed  and 
defined  in  the  Notation  (Appendix  C). 
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In  general,  n  is  greater  than  about  4,  and  figure  2A  indicates  that  (2) 
is  an  adequate  representation  of  the  profiles  for  an  integral  analysis, 
in  which  small  deviations  between  the  actual  and  formulated  profiles 
have  relatively  small  effects.  The  equation  expressing  conservation  of 
moment-of -momentum  (torque)  about  the  centroid  of  the  control  volume 
shown  in  figure  2  is 

d  r  2  .  ,  d 

p  /  J°  7-  (y  -  7)  dr  dy  rd<)>  -  Wp  ^  [  /  uv  (y  -  -|)  dy  d$ 

-Tor!<ro2-r,2>£=°  (3) 


in  which  TQr  =  radial  component  of  the  bed  shear  stress.  The  radial 
shear  force  exerted  on  the  bed  results  principally  from  the  velocity 
profile  just  above  the  bed  being  skewed  by  the  secondary-flow 
velocity.  The  secondary  velocity  itself  is  relatively  small  compared  to 
the  primary  velocity,  and  alone  produces  a  minor  increase  in  the  total 
shear  stress.  It  appears  reasonable  to  assume,  therefore,  that  it  is 
the  skewing  of  the  primary-flow  bed  shear  stress  that  produces  the 
radial  component  of  the  bed  shear,  and  that  the  latter  is  proportional 
to  the  skewing  of  the  velocity  profile.  This  will  be  expressed  as 


*  8  — 
Tos  V 


(4) 


where  0  =  proportionality  factor  to  be  determined  from  measured  rates  of 
streamwise  development  of  secondary  flow  in  curved  channels,  and  Y  = 
mean  streamwise  flow  velocity.  The  quantity  tqs  (hereinafter  denoted 
simply  as  tq)  will  be  related  to  the  local  mean  velocity  by  means  of  the 
Darcy -We is bach  equation. 


where  f  =  Darcy -We Is bach  friction  factor.  Substitution  of  (1),  (2), 
(4),  and  (5)  into  (3)  and  carrying  out  the  indicated  integrations  yields 


dc  dU  .  dc  tt 

in  which 


>c  -  1  <r1  +  ro) 

_  (3n+l)(2n+l)  Bf 
2nSn+l  ° 


g  -  (3n+l) (2n+l) (n+1)  (9) 

2  n(n+2) (2n2+n+l) 

Equation  6  is  a  linear,  ordinary  differential  equation  which  has  for  its 
solution 


U(s)  -  U(S0)  ♦  g2V  exp  C-g1(-jra-)3  /  p-fcy  exp  [g,  (-j-^MlO) 

s0  c 

where  the  change  of  variable  ds  *  Rc  d^  has  been  made.  Note  that  the 
subscript  c  is  used  hereinafter  to  refer  to  centerline  values.  The 
quantity  U(sQ)  is  the  secondary -current  strength  at  s  =  s0.  In  a  field 
application,  the  centerline  curvature,  1/Rc(s),  would  be  determined  from 
a  map  or  survey  and,  in  the  case  of  complex  channel  lineament,  the 
quadrature  appearing  in  (10)  likely  would  have  to  be  evaluated 
numerically.  This  poses  no  problem  inasmuch  as  the  governing  equation 
themselves  must  be  treated  numerically. 


12.  To  determine  the  bed  topography,  and  therefrom  the  streamwise 
and  transverse  distributions  of  flow  depth  for  utilization  in  the 
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numerical  solution  of  the  equations  of  continuity  and  motion  of  the 
fluid,  two  assumptions  will  be  made,  as  follows: 

1.  The  transverse  bed  profile  is  a  straight  line  at  every  channel 
section.  Figures  4A  and  6A,  and  also  the  results  presented  by 
Zimmermann  and  Kennedy®,  demonstrate  that  the  deviations  of 
both  measured  and  more  accurately  computed  transverse  profiles 
from  a  straight  line  are  relatively  small. 

ii.  The  transverse  bed  slope  varies  linearly  with  U.  This  is 
consistent  with  (4)  and  the  bed-layer  equilibrium  analysis 
presented  in  Appendix  A,  where  it  is  shown  that  the  local  bed 
slope  varies  linearly  with  the  local  stress  (14,  App.  A).  In 
terms  of  the  mean  transverse  bed  slope,  Sy(=  sin  b  in  (14, 
App.  A)),  this  equation  reads 

Tor  =  V1_P)  Ap  9  ST  U1) 

where  p  =  bed-layer  porosity;  a p  =  ps_p»  in  which 
p$  *  particle  density  and  p  =  fluid  density;  and 

g  =  gravitational  constant.  Substitution  of  (4),  (5),  and 
(10),  (15,  App.  A),  (16,  App.  A),  and  (17,  App.  A),  into  (11) 
yields 

sT=.Thr'F  -==u  (12> 

f  °50 

where  e  =  Shields  parameter  defined  by  (16,  App.  A)  and  a  = 
proportionality  constant  between  the  bed-layer  thickness,  yj, 
(see  figure  1),  and  shear  velocity,  u*,  used  by  Karim7. 
Equation  (12)  can  be  simplified  to  the  following  expression: 


(12A) 


C  u 

ST  =  93  y 

where 

93  g  4~~:=:7  (12B) 

'9  f  Dso 

Note  that  n  and  f  are  related  by  Nunner's  relation  (17,  App. 

A),  which  is 

n  =  1/vT  (13) 

13.  It  will  be  assumed  further  that  the  depth  changes  across  any 
section  due  to  curvature-induced  inclination  of  the  water  surface  are 
very  small  compared  to  those  due  to  bed  warping.  (Note,  however,  that 
the  effect  of  radial  water-surface  inclination  is  retained  initially  in 
the  equations  of  motion  developed  in  the  next  section,  but  is  then  shown 
to  be  negligible  in  the  streamwise  momentum  equation  for  most  meandering 
river  situations.)  The  depth  at  any  point  across  the  channel  is  then 
given  by 

d(r,s)  *  dc  ±  Syr  (14) 

where  the  sign,  +,  is  adopted  according  as  Rc  (see  figure  3)  is  positive 
or  negative.  The  bed  elevation  at  any  point  in  the  channel  then  is 
given  by 

b(r,s)  =  hc(0,sQ)  -  Sc(s-sQ)  ±  STr  (15) 

in  which,  again,  the  sign  is  chosen  to  be  the  same  as  that  of  Rc» 
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14.  The  steady -flow,  depth-integrated  conservation  equations  for 
mass  and  for  radial  and  streamwise  momenta,  expressed  in  radial 
coordinates,  will  be  used  for  calculation  of  the  velocity  field.  The 
continuity  equation  is 

Is*  [V ( H— h ) ]  +  Ip  [r  TT(H-h)]  =  0  (16) 

in  which  TI  =  shift  velocity  (see  figure  1)  which  accounts  for  the 
transverse  mass  shift  that  occurs  in  channels  with  nonuniform  curvature 
(e.g.,  along  meandering  channels  as  the  thalweg  moves  from  the  vicinity 
of  one  bank  to  the  other). 

15.  The  radial -momentum  equation  is 

H  2 

[/  pp-  dr  dy]  rd$  +  j  pg(H-h)2rd$]dr 

+  pg(H-h)2dr  d$  -  pg(H-h)  dr  rd<j>  +  tor  rd<j»  dr  = 

H  H 

P(u+U)v  dr  dy]d$  +  -fjr  p(u+U)2rd<fr  dy]dr  (17) 

Substitution  of  (1)  and  (2)  into  (17)  yields  for  the  depth-integrated 
radial -momentum  equation 

-  9<H-h>  f  * 

Ij.  [V(H-h)(T  *  2^j)]  +  -i  Ip  [r(H-h)(-§i  +  u2)]  (18) 

16.  The  corresponding  equation  expressing  conservation  of 
streamwise  momentum  is 

tt  Ij  pg(H-h)2dr]d$  -  pg(H-h)  rd$  dr  j  +  T  rd$  dr  = 


(19) 


H  H 

pv2dr  dy]d$  +  -|jr  pv(u+TJ)rd<fr  dy]dr 

which,  after  introduction  of  the  velocity  distributions  adopted  for  this 
analysis,  (1),  (2),  and  the  uniformly  distributed  transverse  shift 

velocity,  yields  the  following  depth-integrated  streamwise-momentum 
equation: 


The  numerical  treatment  of  these  equations  is  described  in  PART  III. 


Sediment-Discharge  Relation 


17.  The  local  sediment  discharge  will  be  calcualted  on  the  basis 
of  the  local  streamwise  velocity  using  a  power-law  relation, 

qt  -  a  Vb  (21) 

in  which  qt  =  total  sediment  discharge  per  unit  width;  and  the 

coefficient  a  and  exponent  b  are  to  be  determined  on  the  basis  of  a 
sediment -discharge  predictor  or  data  on  the  channel  under  consideration 
for  its  particular  flow  regime,  bed-material  size,  etc.,  etc.  The 
numerical  program  is  structured  such  that  other  sediment-discharge 
relations  can  be  incorporated  into  it.  In  particular,  it  is  envisioned 

that  the  future  development  might  utilize  a  formulation  such  as  that 

recently  developed  by  Karim?,  which  uses  an  iterative  procedure  to 
calculate  sediment  discharge  and  friction  factor  as  interdependent 
variables.  This  would  permit  incorporation  of  laterally  nonuniform 
friction  factor  into  the  program,  and  calculation  of  the  sediment 

discharge  of  each  bed-material  size  fraction.  However,  time  and  funds 
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did  not  permit  undertaking  of  this  rather  major  effort  in  the  present 
study. 

18.  It  is  recognized  that  the  nonlinearity  of  (21)  can  lead  to 
calculated  streamwise  variations  in  the  sediment  discharge  along  a 
nonuni formly  curved  channel  with  warped  bed.  A  correction  procedure  is 
incorporated  into  the  numerical  model  which  compensates  for  this 
artifact  in  the  following  way: 

i.  The  sediment  discharges  computed  from  (21)  for  radial 
computational  increments  across  each  computation  section  are 
summed  to  obtain  the  computed  total  sediment  discharge  for  the 
section. 

ii.  The  sediment  discharge  in  each  radial  computational  increment 
is  corrected  by  a  factor  equal  to  the  ratio  of  the  sediment 
discharge  into  the  bend  divided  by  the  computed  total  sediment 
discharge  across  the  section. 

This  insures  that  sediment  continuity  is  preserved  along  the  channel 


PART  III:  NUMERICAL  MODEL 


Numerical  Strategy 

19.  The  three  governing  equations,  (16),  (18),  and  (20),  contain 
three  unknowns:  the  depth-integrated  streamwise  velocity,  V(r,s);  the 
shift  velocity,  TT(r,s);  and  water-surface  elevation,  H(r,s).  The 
secondary-flow  velocity  U(s),  was  calculated  from  the  torsion-balance 
analysis  and  is  given  by  (10),  and  the  bed  elevation,  h(r,s),  was 
obtained  from  the  computed  average  radial  bed  slope  and  expressed  by 
(15).  Numerical  solution  of  these  three  strongly  coupled  equations 
proved  to  be  quite  difficult,  but  was  greatly  simplified  by  introducing 
the  following  restriction.  Note  that  in  (16)  and  (20),  H  appears  only 
in  the  combination  (H-h),  which  is  the  local  depth  given  by  (14),  except 
in  the  first  term  of  (20).  The  streamwise  water-surface  slope  comprises 
two  parts:  one  due  to  the  friction  slope,  which  is  of  order 

(f  )f  -  0(f  £-)  (22) 

and  a  second  resulting  from  the  centrifugal ly-induced  superelevation  of 
the  water  surface  and  of  order 

(f  >S  -  °flr  177)  <23> 

in  which  L  =  characteristic  length  of  the  curve  (say,  the  half¬ 
wavelength  of  a  meander).  The  second  of  these  can  be  neglected  compared 
to  the  first  if 
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«  f 


(24) 


which  is  satisfied  by  most  natural,  sand-bed,  meandering  streams  flowing 
in  the  ripple-  or  dune-bed  regimes.  If  (24)  is  satisfied,  in  (20) 
may  be  replaced  by 
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which  states  that  the  water-surface  elevation  is  constant  across  all 
sections.  Substitution  of  (25)  into  (20)  yields  an  equation  which, 
together  with  (16)  forms  a  pair  of  simultaneous  equations  for  the  two 
velocities  of  interest,  TT  and  V. 

20.  It  is  convenient  for  numerical  analysis  to  simplify  (16)  and 
(20)  as  much  as  possible.  The  radial  coordinate,  r,  is  first  replaced 

by  Rc  +  r  in  (16)  and  the  expressions  for  d(r,s)  and  Sj(s),  (14)  and 
(12),  are  introduced,  which  yields 


Fl“  +  F2V  +  dcF4  <f  +  f  >  -  0 


where  d  +  ST(s)r 

F1  =  M5)  +  -T-T7 -  W 

F2  =  ^2^3  R-  ”  ^1  d-  (28) 

and  C  C 

F4=l+^-ST(s)  (29) 

The  flow-continuity  equation  (26)  is  normalized  using  the  following 
variables: 

IT1  =  — •  v'  =  — •  R1  =  — •  d*  =  — • 

Y  Y  c  w  »  ac  w  » 

r'  =  J;  and  s'  =  J  (30) 

The  normalized  expression  of  (26)  becomes,  after  dropping  the  prime 
superscripts. 


Fiu  *  f2v  *  dcF4  (f +  $  ■ 0 


(26A) 


Similarly,  substitution  of  (14)  into  (20)  and  nondimensionalization  of 
the  equation  yields 


r  r 


Fluv  +  dcF4  It  (uv>  +  'Zn+T  CF1V  +  dcF4 

.  r(n*l)2  r  ,  f,  „2  .  (n+ir  , 

+  F2 +  s]  v  +  oriy  d( 

R. 


3  V- 


2  w2 

f  il. 

V  4  3S 


where 


1  r  r  C 

77  Scp4  OF 
Fr  c 


-2  V2 

>  =F7 


(31) 


(32) 


21.  It  is  also  advantageous  in  the  numerical  treatment  of  the 
equations  to  avoid  computing  small  derivatives  of  the  dependent 
variables.  Because  the  shift  velocity  is  much  smaller  than  the  depth- 
averaged  streamwise  velocity,  the  term,  aU/ar  in  (31)  will  be  eliminated 
by  the  use  of  (26A),  with  the  result 


W  *  2n*T-l  fr  *  2iv*T  UV 

*  dcF4  [irffcjT  *  Ir  *  + I3  s 


.  R 

1  r  r  C 

77  Scf4  1 TTF 
Fr  c 


(31A) 


22.  The  numerical  strategy  employed  to  solve  for  TJ  and  V 
proceeded  as  follows: 

i.  The  local  depth-averaged  velocity  was  approximated  by  the 
Darcy-Weisbach  equation. 


(33) 


in  which  d  is  given  by  (14)  and  S  by 


s  =  s 


Rc 

c  FTF 
c 


which  expresses  continuity  of  energy  slope  across  the  channel. 

ii.  The  velocity  V  given  by  (33)  was  substituted  into  (16)  which 
was  then  integrated  numerically  to  obtain  the  first  estimate 
for  IT. 

iii.  The  value  of  TT  was  substituted  into  ( 3 1A )  which  was  integrated 
to  obtain  the  next  estimate  for  V. 

iv.  The  V  computed  in  step  iii  was  substituted  into  (16)  and  a  new 
estimate  of  TJ  obtained. 

v.  The  iteration  procedure  between  (16)  and  (31A)  was  continued 
until  satisfactory  convergence,  as  measured  by  the  differences 
between  successive  values  of  TT  and  V,  was  obtained. 

Further  details  are  presented  below. 


Numerical  Solutions  for  TT 

23.  In  order  to  solve  (26A)  and  (31A)  numerically  for  the  two 
unknown  variables,  TI  and  V,  a  backward  finite-difference  scheme  was 
employed.  Figure  4  shows  the  coordinate-grid  layout  that  was  utilized. 
The  indicies  i  and  j  represent  the  streamwise  and  radial  positions, 
respectively.  Note  that  the  origin  of  the  transverse  coordinate  was 
taken  at  the  channel  centerline.  In  discretizing  both  radial  and 
strearrwise  derivatives  of  an  arbitrary  variable  F,  the  following 
backward  finite-difference  scheme  was  utilized: 

1L  .  FVj  Fi.J-1  (35) 

3r  rj  rj-l 
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and 


3F  _  Fi,j  "  Fi -l,j 
as  "  s.  .  -  s.  ,  . 
i.j  i-i,  j 


24.  An  approximate  solution  for  U  can  be  obtained  from  the  flow- 
continuity  equation,  (16),  by  introducing  the  Darcy-Wei sbach 
relationship  and  (34): 


,2  89dSc  Rc 


Substitution  of  (37)  into  (16),  use  of  (14)  and  (12A),  and  subsequent 

discretization  yields  the  following  explicit  expression  for  TJ. : 

J 


_  _  «<V)|,.i 

uj  ■  vi  dp^rir 


srsi  j 

2isJ  exp["9i  ~Tc  ] 

3g2g3  /  1  d(Rr+r)| •  T 

^  J 


where 


2  ^c3T_<*c  ^1 

T  =  2Rc2(RcSt  -  dc)[  i  -c- 

(5RcVdc>T2  -  (dc  <■  3RcSt)T3. 


T>  ' 


T  4-r  r\  oT 

h  =  C  T 


T?  =  7IrT  'j-1 

t-/RT  c  T 

- - -  in  I - —I :  R  ST  >  0 

2/irsr  't+/rsr'  c  ' 


— \ - tan"1  [ - - — ]:  R  ST  <  0 

/-rcst  /-rcst 
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t  =  / 


Rc(dc  +  Syr) 


Note  that  the  boundary  condition  TI  =  0  was  imposed  at  the  inside 
(convex)  bank.  Equation  38  gives  the  approximate  solution  for  U  which 
can  be  substituted  into  the  streamwi se -momentum  equation,  (31A),  to 
solve  for  V. 

25.  Once  V  is  computed  at  Section  I  =  i,  the  flow-continuity 
equation  can  be  again  utilized  to  solve  for  TJ.  .  in  the  iterative 
process  without  utilizing  the  Darcy -Wei sbach  relationship.  The  final 
discretized  form  of  (16)  in  terms  of  V  is 

*0  -  Vi  *1*73 


c(Rc+d)^|  t  - 


<  +r)  . 
c 


Numerical  Solution  for  V 

26.  Discretization  of  the  streamwise-momentum  equation,  (31A), 
yields  the  following  quadratic  equation  for  : 


where 


AV^  •  +  BV i  .  +  C  =  0 


F,  -  d  F.  -i  -t 

“  n(n+2)  8  Sisj  -  si_1j  Ln(n  +  2)  2J 

d  F  U.  .  F 

8  -  ft-M- t  [U1 , j  +Sii]  *  2HTTuu 


(48) 


dCF4  r  1  ,  v2  J_  c  r  RC 
-  Sisj  -  Sl„ltj  l^2j  2^  Vi-l,j  -  p2  bcF4^TF- 

It  should  be  noted  that  the  total  water  discharge  calculated  with  the 
computed  transverse  distribution  of  V  did  not  equal  the  imposed  total 
discharge  due  to  discretization  errors.  Therefore,  an  adjustment  was 
made  to  V  at  each  cross  section  by  multiplying  j  by  the  ratio  of  the 
imposed  water  discharge  to  the  computed  water  discharge.  This  ratio  was 
typically  of  the  order  of  1.0005. 


Boundary  Conditions 

27.  The  streamwise  velocity,  V,  was  specified  at  the  inlet 
section  (s  =  sQ)  and  along  the  inside  bank  of  the  computational  reach  by 
the  Darcy -Wei sbach  relationship,  (37).  Along  the  inside  bank,  the  shift 
velocity,  XT,  was  set  equal  to  zero. 


Computer  Program 

28.  The  program  PR-SEG6  consists  of  a  main  program,  four 
subroutines,  and  seven  sub-subroutines.  Listings  of  the  main  program, 
the  subroutines,  a  sample  input  file,  and  a  sample  output  file  are 
included  in  Appendix  B.  Note  that  the  sample  input  and  output  shown  in 
Appendix  B  are  for  the  idealized  two-bend  model  which  is  discussed  in 
PART  IV. 

29.  The  main  program  first  reads  the  common  input  variables  from 
the  input  file  called  SEGDAT:  V,  dc,  W,  Sc.  P»  p  ,  and  NSEG  (number  of 
channel  segments  in  the  reach  that  requires  new  input  parameters).  The 
program,  then  reads  the  following  parameters  at  each  new  channel 
segment:  a,  b,  M  (number  of  radial  positions),  N  (number  of  streamwise 
positions),  Rc,  s0,  si  (centerline  streamwise  coordinate  of  the 
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downstream  end  of  the  segment),  a,  0,  9c>  D^,  and  NSTEPS  (number  of 
cross  sections  into  which  the  channel  segment  is  divided).  The  program 
computes  the  boundary  values  of  U  and  V  across  the  inlet  section  using 
(10)  and  (37),  and  the  transverse  distribution  of  U  is  subsequently 
computed  from  (38).  The  program  then  advances  to  the  next  downstream 
section,  and  computes  U  and  V  according  to  the  iterative  scheme 
described  in  paragraphs  22,  24,  25,  and  26. 

30.  Subroutine  PQN  was  used  to  determine  the  total  water 
discharge  for  a  given  cross  section  with  computed  transverse 
distributions  of  streamwise  velocity  and  depth.  Subroutine  PG 
determines  the  g-parameters  defined  by  (8),  (9),  and  (12B).  Subroutine 
EVAL  evaluates  the  shift  velocity,  TT,  given  by  (38). 

31.  The  main  program  writes  the  following  outputs  in  the  output 
file,  called  OUTT,  for  each  cross  section:  Sy,  uc,  number  of  iterations 
required  to  compute  satisfactory  convergence  of  "u  and  V,  and  transverse 
distributions  of  IT,  V,  d,  U  +  H,  and  qt. 


Sensitivity  Analysis 


32.  The  effects  of  the  specified  error  tolerance  for  IT,  the  grid 
size,  the  transverse  derivative  of  U,  and  parameters  «  and  0  on  TT  and  V 
were  tested  using  the  basic  hydraulic  and  sediment  parameters  that  were 
utilized  in  the  Oakdale  flume  experiments  conducted  at  the  Iowa 
Institute  of  Hydraulic  Research,  The  University  of  Iowa,  by  Odgaard  and 
Kennedy3.  The  basic  parameters  were:  7  =1.56  ft/s,  dc  =  0.505  ft,  W  = 
8.0  ft,  Rc  =  43  ft  (see  figure  5),  Sc  =  0.00104,  p  =  0.4,  Ps/P  = 
2.65,  v  =  1.21  x  10"5  ft^/s,  D50  =  0.3  mm,  and  9C  =  0.032. 

33.  The  relative  errors  for  TT  and  V  computed  at  each  cross 
section  were  defined  by 
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where  aTJ  and  aV  are  changes  in  U  and  V  between  adjacent  radial 
positions,  respectively.  Because  TJ  was  typically  at  least  two  orders  of 
magntiude  smaller  than  V,  the  error  tolerance  for  TJ  was  selected  to  be 
an  order  of  magnitude  larger  than  that  for  V.  In  the  sensitivity  tests, 
parameters  a  and  8  were  set  equal  to  1.0  and  3.35,  respectively,  and  a 
grid  size  of  6  in.  was  used.  For  Ey  of  0.1%,  two  tests  were  run 
for  E-g  values  of  2%  and  0.4%.  There  were  no  discernible  differences 
between  two  sets  of  values  computed.  An  additional  run  with  Eg  equal  to 
0.2%  demonstrated  that  this  criterion  was  not  able  to  be  satisfied  with 
single-precision  computations.  Note  that  Ey  was  of  the  order  of  10“5 
between  successive  iterations  for  V.  It  was  concluded  that  satisfactory 
results  could  be  obtained  with  the  error  tolerances  of  Eg  and  Ey  equal 
to  approximately  2%  and  0.1%,  respectively. 

34.  Sensitivity  tests  were  run  for  different  square-grid  sizes. 
The  grid  size  was  reduced  step  by  step  until  no  significant  changes  in 
estimates  of  U  and  V  resulted.  As  shown  in  figures  6  and  7,  grid  sizes 
of  4  in.  and  6  in.  yielded  quite  similar  transverse  distributions 
of  U  and  V;  in  fact,  the  two  sets  were  almost  identical  at  the 
downstream  end  of  the  channel  bend.  From  the  sensitivity  analysis,  it 
was  concluded  that  the  grid  size  should  be  approximately  equal  to  the 
mean  flow  depth.  Note  that  the  mean  flow  depth  of  the  Oakdale  flume  was 
about  6  in.  (see  figure  5). 

35.  In  obtaining  the  simplified  streamwise  momentum  equation, 
(31),  the  secondary -flow  velocity,  U,  was  treated  as  a  function  of  only 
s  because  of  the  assumption  of  constant  transverse  bed  slope,  as  given 
by  (12).  However,  in  computing  V  by  means  of  (45),  the  computer  program 
utilized  a  radi al ly-varied  secondary-flow  velocity  distribution  derived 
by  Falcon  and  Kennedy^ 


U  -  /V  wd  w  c  v 
IT  "  (V.)(d  ^F7r' 


(50) 
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was  utilized.  When  (50)  is  substituted  into  (20),  the  discretized 
st reamwise -momentum  equation,  (45),  yields  coefficients  A,  B,  and  C  that 
are  slightly  different  from  those  given  by  (46),  (47),  and  (48).  In 
order  to  ascertain  the  validity  of  the  computational  scheme,  a  special 
test  run  was  made  with  the  term,  3U/3r,  retained  in  the  streamwise 
momentum  equation.  The  computed  distributions  of  U  and  V  are  compared 
in  figures  8  and  9  with  those  obtained  using  (45),  which  was  developed 
without  the  term,  3U/3r.  As  can  be  seen  in  these  figures,  the  effects 
of  the  term,  3U/3r,  on  overall  estimates  of  IT  and  V  are  minor. 

36.  The  parameters  a  and  3  control  the  transverse  bed  slope,  Sy» 
and  the  development  rate  of  the  secondary-flow  velocity,  respectively, 
as  can  be  seen  from  (12),  and  (8)  and  (10).  Figures  10  and  11  show  the 
effects  of  a  on  the  transverse  distributions  of  TT  and  V.  As  can  be  seen 
in  these  figures,  the  smaller  a  resulted  in  larger  Sj,  and  consequently 
in  much  larger  shift  velocities  along  the  initial  entrance  reach  of  the 
bend.  The  smaller  a  also  resulted  in  much  smaller  streamwise  velocities 
along  the  inside  bank,  because  the  larger  Sy  decreased  the  flow  depth 
there.  Similar  effects  of  3  on  U  and  V  are  seen  in  figures  12  and  13. 
The  smaller  3  resulted  in  a  slower  development  rate  of  the  secondary- 
flow  velocity,  and  reduced  the  rate  of  the  development  of  the  transverse 
nonuniformity  in  V. 


PART  IV:  RESULTS  OF  NUMERICAL  SIMULATIONS 


Oakdale  Flume 

37.  The  Oakdale  flume  shown  in  figure  5  is  a  l:48-scale,  highly 
idealized,  undistorted  model  of  the  Sacramento  River  bend  lying  between 
R.M.  188  and  189,  approximately.  Experimental  data  on  the  streamwise 
distribution  of  the  equilibrium  transverse  bed  slope  and  transverse 
distributions  of  the  depth-averaged  streamwise  velocity  reported  by 
Odgaard  and  Kennedy^  were  compared  with  the  computer-simulated 
results.  The  basic  hydraulic  and  sediment  parameters  described  in 
paragraph  32  were  utilized  in  the  simulation.  Additional  parameters 
specified  were:  a  =  1.00,  a  =  3.35,  grid  size  =  6  in.,  n  =  4.24,  E-g  = 
2%,  and  Ev  =  0.1%. 

38.  Figures  14  and  15  demonstrate  generally  good  agreements 
between  the  measured  and  computed  transverse  distributions  of  the  flow 
depth  and  the  depth-averaged  streamwise  velocity,  respectively,  at  ,|>  = 
20°.  The  results  for  ^  =  114°  are  shown  in  figures  16  and  17,  in  which 
the  observed  streamwise  velocities  are  seen  to  be  somewhat  larger  than 
the  computed  values  in  the  outside  portion  of  the  channel.  The  larger 
measured  velocities  near  the  outside  bank  are  believed  to  be 
attributable  to  the  very  low  roughness  of  the  exposed  plywood  bank  of 
the  trapezoidal  flume  section.  Note  that  the  friction  factor  was  kept 
constant  in  the  whole  flow  field  in  the  numerical  model.  Figure  18 
depicts  extremely  good  agreement  between  the  measured  and  computed 
transverse  bed  profiles  for  $  =  146°. 

Sacramento  River 

39.  The  Sacramento  River  bend  between  R.M.  188  and  189  shown  in 
figure  19  was  simulated  for  two  water  discharges  (Q  =  9,000  cfs  and 
25,800  cfs).  Basic  field  data  were  collected  in  the  reach  in  1979  and 


1980  by  the  U.S.  Geological  Survey  (USGS)  (Odgaard  and  Kennedy3).  As 
can  be  seen  in  table  1,  both  the  hydraulic  and  sediment  parameters 
varied  widely  along  the  bend.  Therefore,  average  values  of  the  various 
quantities  listed  in  table  1  were  utilized  for  the  numerical 
simulations. 

40.  The  transverse  distributions  of  the  measured  and  computed 
flow  depth  and  streamwise  velocity  for  Q  =  9,000  cfs  are  shown  for  $  = 
80°  and  126°  in  figures  20  and  21,  and  figures  22  and  23, 
respecti vely.  Note  that  at  <j>  =  80°,  velocities  and  depths  were  measured 
at  only  four  verticals  across  the  channel,  while  they  were  measured  at 
ten  verticals  at  <j>  =  126°.  Despite  the  fact  that  averaged  input  data 
were  adopted  for  the  simulation,  the  numerical  model  reproduced  the 
field  distributions  surprisingly  well  for  the  low  river  discharge  of 
9,000  cfs.  The  distributions  obtained  for  the  higher  discharge  of 
25,800  cfs  are  shown  in  figures  24  through  27.  The  agreements  between 
the  measured  and  predicted  values  are  seen  to  be  not  as  good  as  those 
for  Q  =  9,000  cfs;  however,  it  is  believed  that  during  high  flows  the 
channel  bed  had  not  attained  an  equilibrium  configuration.  For  example, 
the  measured  transverse  bed  slopes  shown  in  figures  22  for  Q  =  9,000  cfs 
and  figure  26  for  Q  =  25,800  cfs  are  entirely  different.  The  field 
transverse  bed  slope  was,  paradoxically,  much  smaller  during  the  high 
flow,  resulting  in  the  decreasing  streamwise  velocity  toward  the  outside 
bank,  as  seen  in  figure  27.  This  type  of  abnormal  transient  phenomenon 
likely  is  a  consequence  of  the  rapidly  changing  flow  conditions,  and 
cannot  be  simulated  by  a  steady-state  numerical  model.  It  should  be 
noted  that  each  Sacramento  River  simulation  required  approximately  0.7 
second  CPU  time  per  100-grid  points  using  the  PRIME-750  computer  at  The 
University  of  Iowa. 


41.  The  numerical  results  presented  in  paragraphs  37  through  40 
for  the  Oakdale  flume  and  the  Sacramento  River  were  obtained  using 
constant  centerline  curvature.  In  order  to  demonstrate  the  ability  of 
the  computer  program  to  handle  nonuniform  curvature,  two  simulations 
were  made  for  single  bends  with  gradually  varying  centerline 
curvature.  These  numerical  simulations  were  made  also  to  illustrate  the 
behavior  of  flow  in  idealized,  nonuniform  river  bends. 

42.  The  first  simulation  was  for  a  four-segment  channel  bend  with 
stepped  decreases  in  curvature  in  the  downstream  direction,  as  depicted 
in  figure  28.  The  centerline  radius  of  the  first  segment  was  2,000  ft, 
and  this  value  was  increased  by  2.5%  for  each  of  the  subsequent  three 
segments,  resulting  in  a  total  channel  length  of  7,400  ft.  It  was  found 
that  a  5%  increase  in  Rc  produced  such  large  transverse-bed-slope 
changes,  which  appear  as  sloped  steps  in  the  bed  elevation,  that  the 
program  would  not  run.  Therefore,  in  cases  in  which  Rc  increases  along 
a  bend,  the  curve  should  be  subdivided  into  sufficiently  short 
subreaches  that  the  increments  in  Rc  are  less  than  about  2.5%,  although, 
as  discussed  in  the  next  example,  the  model  can  accommodate  larger 
changes  in  the  case  of  decreasing  Rc,  The  basic  hydraulic  and  sediment 
parameters  used  were  identical  to  those  for  the  Sacramento  River  at  high 
flow,  listed  in  table  1.  A  grid  size  of  14.5  ft  was  used,  and  the 
parameters  a  and  3  were  set  at  0.86  and  7.13,  respectively.  The 
computed  longitudinal  and  transverse  distributions  of  the  normalized 
shift  velocity  are  shown  in  figures  29  and  30,  respectively.  In  figure 
29,  the  shift  velocities  computed  for  sections  65,  193,  321,  and  449  are 
connected  by  straight  lines.  The  shift  velocity  developed  rapidly  in 
the  first  segment,  with  its  maximum  values  occurring  near  r/W  equal  to 
-0.25,  and  diminished  gradually  after  section  193.  At  section  385,  the 
shift  velocity  along  r/W  equal  to  -0.25  became  negative,  and  remained  so 
until  section  469.  This  flow  redistribution  directed  radially  inward 
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was  a  consequence  of  the  increased  Rc.  Figures  31  and  32  depict  the 
longitudinal  and  transverse  distributions  of  the  depth-averaged 
streamwise  velocity,  respectively.  Along  the  inside  bank,  the 
streamwise  velocity  decreased  initially;  however,  it  increased  farther 
downstream  as  the  larger  radii  of  curvature  produced  less  steep 
transverse  bed  slopes.  The  values  of  Sy  at  sections  1,  65,  193,  321, 
449,  and  513  were  0,  0.058,  0.063,  0.062,  and  0.060,  respectively. 

43.  The  second  idealized  case  simulated  was  a  single  bend  with 
radius  of  curvature  that  decreased  10%  between  curve  subreaches.  The 
numerical  results  are  not  presented  herein,  because  the  qualitative 
characteristics  are  very  similar  to  those  for  the  two-bend  curve  with 
decreasing  radius  of  curvature  presented  in  the  following  section. 
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V. 


Idealized  Two-Bend  Model  with  Gradually 
Decreasing  Radius  of  Curvature 


44.  An  idealized  two-bend  model,  shown  in  figure  33,  was 
tested.  The  two-bend  reach  consisted  of  four  segments  with  equal 
centerline  length  of  67.5  ft.  The  centerline  radius  of  curvature  of  the 
first  segment  was  43.0  ft,  and  was  reduced  by  10%  for  each  subsequent 
subreach.  The  sign  of  Rc  was  reversed  after  the  second  subreach.  The 
simulation  was  made  on  the  basis  of  the  principal  parameters  used  in  the 
Oakdale  flume  simulation.  These  parameters  are  described  in  paragraph 
37,  except  that  a  =  1.42  and  0  =  3.28  were  used  in  the  present 
simulation.  Figures  34  and  35  show  the  longitudinal  and  transverse 
distributions  of  the  normalized  shift  velocity,  TJ/V,  respectively.  The 
shift  velocity  increased  rapidly  in  the  first  segment  and  decreased 
gradually  toward  the  end  of  the  first  bend.  Once  the  flow  entered  the 
second  bend,  a  mass  shift  took  place  toward  the  right  bank  due  to  the 
change  in  sign  of  the  channel  curvature.  Note  that  in  figure  34,  the 
computed  data  points  at  sections  69,  205,  341,  477,  and  545  are 

connected  by  straight  lines.  As  shown  in  figure  35,  the  maximum  value 
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of  the  shift  velocity  across  the  cross  section  was  closer  to  the  convex 
side  of  the  bend.  Figures  36  show  the  transverse  distributions  of  the 
depth-averaged  streamwise  velocity  computed  at  sections  1,  205,  341,  and 
545.  The  transverse  location  of  the  maximum  V  gradually  shifted 
radially  outward  in  the  first  segment,  and  reached  the  outside  bank  at 
section  71.  The  maximum  V  remained  along  the  left  bank  until  section 
273,  after  which  the  flow  became  concentrated  near  the  right  bank.  The 
maximum  V  reached  the  right  bank  at  section  417  in  the  second  bend. 
Because  the  streamwise  velocity  at  the  left  bank  at  section  273  was  much 
larger  than  that  at  section  1,  a  larger  streamwise  distance  was  required 
to  attain  redistribution  of  the  flow  in  the  second  bend. 

45.  Figure  37  shows  the  transverse  distributions  of  the  unit 
total-load  discharge,  q^,  computed  at  various  cross  sections.  The 
sediment -transport  coefficients  a  and  b  in  (21)  were  taken  to  be  0.108 
and  4.0,  respectively.  Note  that  the  units  of  V  and  qt  are  ft/s  and 
tons/ft/day,  respectively.  These  coefficients  yielded  a  mean  total -load 
concentration  of  300  mg/1  (or  about  5  tons/day)  for  the  Oakdale  flume. 
The  distribution  curves  shown  in  this  figure  are  seen  to  be  generally 
congruent  with  those  transverse  distributions  of  the  streamwise  velocity 
shown  in  figure  36,  because  of  the  sediment -transport  relation  adopted 
being  a  power  function  of  V. 


PART  V:  SUMMARY  AND  CONCLUSIONS 


Conclusions 

46.  The  principal  features  of  the  numerical  model  developed 
herein  for  calculation  of  flow  and  sediment-transport  distributions  in 
alluvial -river  bends  may  be  summarized  as  follows: 

i.  The  secondary -flow  strength  and  the  bed  topography  are 
uncoupled  from  the  calculation  of  distributions  of  lateral 
shift  velocity  and  streamwise  velocity.  This  is  accomplished 
by,  first,  calculating  the  secondary-flow  strength  on  the 
basis  of  conservation  of  flux  of  moment  -of -momentum,  and, 
second,  determining  the  bed  topography  on  the  basis  of  radial 
force  equilibrium  of  the  moving  bed  layer. 

ii.  The  distributions  of  lateral  shift  velocity  and  depth-averaged 
streamwise  velocity  are  calculated,  for  the  warped  channel 
determined  as  described  in  step  i  above,  from  the  depth- 
integrated  equations  expressing  conservation  of  mass  and 
momentum.  It  was  concluded  that  for  flows  which  satisfy  (24), 
it  is  not  necessary  to  include  the  third  conservation 
equation,  that  for  radial -di recti  on  momentum,  or  to  iterate 
among  three  equations  to  obtain  a  solution.  The  numerical 
scheme  utilizes  the  backward  finite-difference  method,  and 
evaluates  transverse  and  streanwise  distributions  of  the 
radial  mass-shift  velocity  and  the  depth-averaged  streamwise 
velocity. 

47.  Numerical  simulations  utilizing  the  model  developed  were  made 
for  one  laboratory  flow,  two  Sacramento  River  flows,  and  three  different 
idealized  channel  bends.  The  principal  conclusions  obtained  from  the 
simulations  are  as  follows: 


i.  Generally  satisfactory  agreement  between  computed  and  measured 

results  was  obtained  by  utilizing  error  tolerances 

of  E-g- and  Ey  of  2%  and  0.2%,  respectively.  In  the  absence  of 
better  information,  it  is  recommended  that  a  =  1.00  and  3  - 
3.50  be  utilized.  In  instances  where  actual  field  data  are 
available  on  the  rate  of  development  and  equilibrium  values  of 
Sy»  a  and  0  should  be  adjusted  on  the  basis  of  the  data. 

ii.  The  most  cost-effective  square-grid  size  is  approximately 
equal  to  the  mean  flow  depth. 

iii.  The  computer  program  is  capable  of  simulating  flow  in 

multiple-bend  channels  with  stepwise-varying  radius  of 

curvature.  On  the  basis  of  the  numerical  simulations,  it  was 
found  that  the  maximum  permissible  stepwise  change  of 
centerline  curvature  for  which  the  program  will  run  is  about 
2.5%  in  the  case  of  increasing  Rc,  and  about  10%  for 
decreasing  Rc. 

48.  Further  development  and  improvement  of  the  model  should 
include  the  following: 

i.  More  complete  and  modern  sediment -di scharge  and  friction- 
factor  models  should  be  incorporated  into  the  model.  In 
particular,  it  is  recommended  that  Karim's?  model  be 
incorporated  into  the  program  to  permit  calculation  of  lateral 
and  streamwise  variations  of  friction  factor  based  on  local 
flow  depth,  velocity,  and  sediment  discharge.  Karim's  model 
is  unique  in  that  it  formally  takes  into  account  the 
interdependence  between  sediment  discharge  and  friction 
factor,  an  interdependency  which  appears  to  be  very  important 
in  channel -bend  flows. 
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ii.  A  further  refinement  of  the  flow  calculation  would  involve 

incorporation  of  the  radial -momentum  equation,  (18),  This 

would  permit  application  of  the  model  to  bends  with  relatively 
short  radius  of  curvature.  However,  the  numerical  model  would 
become  much  more  complex,  and  would  require  significantly  more 
computer  time.  The  model  developed  herein  is  believed  to  be 
adequate  in  its  flow-calculation  aspects  for  most  natural 
al luvial -channel  bends. 

iii.  An  effort  should  be  made  to  incorporate  features  into  the 

model  to  permit  prediction  of  the  occurrence  and 
characteristics  of  point  bars  and  their  effects  on  the  flow 
field.  It  is  believed  that  this  likely  will  require 
incorporation  of  the  radial -momentum  equation  and  a  more 

refined  sediment-discharge  predictor,  as  described  above. 

iv.  As  is  generally  the  case  in  river-flow  analysis,  there  is  a 

pressing  need  for  detailed,  diagnostic-quality  data  on  the 

distributions  of  velocity  and  sediment  discharge  from  both 
natural  and  laboratory  streams. 

v.  After  some  experience  is  gained  with  the  model,  the  computer 

program  should  be  reviewed,  made  more  compact  and  concise 

where  possible,  and  a  user's  manual  for  the  program  should  be 
prepared. 
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Table  1 

Hydraulic  and  Sediment  Parameters  Used  in  Simulating  the  Sacramento  River 


Parameter 


Low  Flow 


High  Flow 


Measured 

Range 


Average 

Value 

Used 


Measured 

Range 


Average 

Value 

Used 
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Figure  8  Transverse  distributions  of  U/V  for  cases  with  and  without  9U/9r  term  in  the 
streamwise  momentum  equation 
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Transverse  distributions  of  measured  and  computed  d/d  for  the  Oakdale  flume 


Figure  18  Transverse  distributions  of  measured  and  computed  d/d  for  the  Oakdale  flume 


Figure  20  Transverse  distributions  of  measured  and  computed  d/d  for  the  Sacramento 
River  at  low  flow  (<p  =  80°)  c 


Transverse  distributions  of  measured  and  computed  V/V  for  the  Sacramen 
River  at  low  flow  (0  =  80°) 


Transverse  distributions  of  measured  and  computed  d/d  for  the  Sacrament 
River  at  high  flow  (<t>  =  80°)  c 


Transverse  distributions  of  measured  and  computed  V/V  for  the  Sacramento 
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Figure  26  Transverse  distributions  of  measured  and  computed  d/d  for  the  Sacramento 
River  at  high  flow  (<}>  =  126°) 


Figure  28  Idealized  single  bend  with  gradually  increasing  radius  of 
centerline  curvature 


grad 


Figure  30  Transverse  distributions  of  computed  U/V  for  ideali 
gradually  increasing  radius  of  centerline  curvature 


Longitudinal  variations  of  computed  V/V  for  idealized  single  bend  with 
gradually  increasing  radius  of  centerline  curvature 


Figure  32  Transverse  distributions  of  computed  V/V  for  idealized  single  bend  with 
gradually  increasing  radius  of  centerline  curvature 


SECTION  I.D.  NUMBER  (I) 

Figure  34  Longitudinal  variations  of  computed  U/V  for  idealized  two-bend  channel 
with  gradually  decreasing  radius  of  centerline  curvature 


Figure  36  Transverse  distributions  of  computed  V/V  for  idealized  two-bend  channel  with 
gradually  decreasing  radius  of  centerline  curvature 


Figure  37  Transverse  distributions  of  computed  q  for  idealized  two-bend  channel 
gradually  decreasins  radius  of  centerline  curvature 
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FLOW  IN  ALLUVIAL-RIVER  CURVES 


by 

Marco  Falcon  Ascanio* 
and 

John  F.  Kennedy*- 

I.  INTRODUCTION 

Even  casual  observers  of  Earth's  geological  features  soon  notice  that 
natural  alluvial  streams  are  seldom  straight  along  reaches  of  more  than  a  few 
channel  widths.  Hydraulic  engineers  and  other  students  of  fluvial  processes 
long  have  recognized  meandering  to  be  not  only  an  intriguing  geometrical  and 
kinematical  feature  of  rivers,  but  also  one  that  has  major  effects  on  their 
sediment-transport  and  roughness  characteristics.  Fluid  mechanicians 
appreciate  further  that  the  internal  structure  of  flow  in  meandering  rivers  is 
as  fascinating  as  their  migrating,  serpentine  channel  lineament.  Especially 
engaging  is  the  interaction  between  the  vertical  profile  of  the  primary  flow 
and  the  centrifugal  forces  resulting  from  the  flow's  curvature.  The  principal 
result  is  the  well  known  spiraling  or  secondary  flow  in  planes  normal  to  the 
channel  axis.  Because  the  bed-surface  sediment  of  a  stream  actively 
transporting  its  bed  material  is  in  a  quasi-fluidized  state,  even  the 
relatively  small  radial  component  of  the  bed  shear  stress  and  small  near-bed 
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velocities  produced  by  the  secondary  flow  transport  sediment  toward  the  inner 
(convex)  banks  until  the  bed  becomes  inclined  such  that  the  gravity  and  shear 
forces  exerted  radially  along  the  bed  on  the  moving  bed-load  particles  are  in 
balance.  The  resulting  greater  depth  near  the  outer  banks  increases  the 
primary-flow  velocities  there,  which  in  turn  intensifies  the  erosive  attack  on 
the  banks  and  also  undermines  them.  Both  of  these  effects  aggravate  bank 
erosion  and  thereby  promote  further  growth  of  the  meanders. 

Although  the  seemingly  disproportionate  effects  of  channel  meandering  on 
river  flow  have  been  appreciated  for  several  decades,  attempts  to  develop  a 
mathematical  model  for  the  secondary  flow  and  its  interactions  with  the 
primary  flow  and  sediment  motion  have  met  with  only  limited  success.  The 
principal  stumbling  block  encountered  arises  from  the  radial  shear-stress 
force  exerted  on  an  elemental  control  volume  at  any  elevation  (the  vertical 
distribution  of  which  is  the  principal  determinant  of  the  radial-plane 
velocity  profile)  being  the  small  difference  between  two  much  larger 
quantities:  the  centrifugal  body  force  and  the  radial  pressure  force 

resulting  from  superelevation  of  the  water  surface.  It  is  important  to  bear 
in  mind  that  even  though  the  integrals  of  these  two  forces  over  the  depth  are 
very  nearly  equal,  locally  they  are  grossly  out  of  balance.  The  radial 
gradient  of  pressure  resulting  from  the  transverse  inclination  of  the  free 
surface  is  very  nearly  constant  over  the  depth,  while  the  centrifugal  force 
varies  from  zero  at  bed  level  to  a  maximum  near  the  free  surface.  In  fact,  it 
is  precisely  the  difference  between  the  distributions  of  these  two  nearly 
equal  forces  that  is  responsible  for  the  secondary  flow.  Moreover,  the 
secondary  flow  (or,  viewed  differently,  the  vertical  gradient  of  the  primary 
velocity)  causes  the  radial  water-surface  slope  to  be  greater  than  it  would  be 


for  a  flow  with  vertically  uniform  primary  velocity  (which  would  not  produce  a 
secondary  current).  This  is  because  the  secondary  flow  produces  an  inward 
radial  shear  force  on  the  bed;  the  corresponding  radially  outward  force  on  the 
flow  must  be  balanced  by  part  of  the  radial  pressure-gradient  force.  Thus, 
just  in  determining  the  distributions  of  the  three  principal  radial  forces-- 
pressure,  shear,  and  centrifugal--exerted  on  the  flow,  one  is  faced  with  the 
problem  of  having  two  of  them--shear  and  pressure--unknown,  even  if  the 
velocity  distribution  of  the  primary  flow  and  hence  also  the  centrifugal-force 
distribution  are  known.  Clearly  to  proceed  with  the  calculation  of  these 
forces,  another  relation,  in  addition  to  the  equation  expressing  the  balance 
of  radial  forces,  is  needed.  Further  physical  considerations  or  assumptions 
must  be  introduced  to  calculate  the  distribution  of  the  radial  velocity. 

In  the  analytical  model  developed  herein  for  vertical  distributions  of 
radial  shear  stress  and  velocity,  and  radial  distributions  of  depth  and 
streamwise  velocity,  an  expression  for  the  conservation  of  moment-of-momentum 
is  the  additional  relation  utilized  to  close  the  formulation  of  the  radial 


forces.  This  aspect  of  the  analysis  is  similar,  for  example,  to  use  of 
equations  expressing  balances  of  forces  and  moments  to  calculate  the 
supporting  forces  on  a  loaded,  simply  supported  beam.  One  of  the  two  unknown 
forces  does  not  appear  in  the  formulation  of  moments  about  one  of  the 
supports,  and  therefore  the  other  can  be  calculated  directly.  A  roughly 
parallel  approach  is  followed  in  the  present  analysis.  Formulation  of  the 
flux  of  moment-of-momentum  about  the  longitudinal  axis  at  the  bed  surface 
yields  an  expression  for  the  radial  pressure  gradient.  The  radial  momentum 
equation  then  is  used  to  obtain  the  vertical  distribution  of  radial  shear 
stress.  The  transverse  bed  profile  is  determined  from  consideration  of  the 
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force  balance  for  the  moving  bed-load  particles.  The  radial-velocity  profile 
is  calculated  by  introducing  into  the  radial  momentum  equation  a  linear 
primary-shear-stress  distribution  and  the  eddy-viscosity  distribution  obtained 
from  the  power-law  distribution  utilized  for  the  primary  velocity.  Finally, 
the  radial  distribution  of  local  depth  corresponding  to  the  bed  profile  is 
used  in  the  calculation  of  the  transverse  distribution  of  depth-averaged 
streamwise  velocity.  The  analysis  is  limited  to  a  channel  of  constant 
centerline  radius,  which  is  a  good  approximation  to  extended  reaches  of  bows 
of  many  strongly  meandering  natural  channels.  Extension  of  the  analytical 
model  to  weakly  meandering  channels  is  developed  in  Falcon's  thesis  (1979). 

Some  of  the  background  literature  on  this  problem  is  cited  in  connection 
with  development  of  the  present  model.  For  a  more  complete  review,  reference 
is  made  to  the  surveys  by  Callander  (1968,  1978),  to  Falcon's  (1979)  thesis, 
and  to  Odgaard's  (1981)  paper  on  river-bend  topography. 

II.  AMALYSIS 

General .  The  idealized  channel  treated  here  has  infinite  length, 
constant  width,  an  erodible  sediment  bed,  and  banks  with  a  common  center  of 
curvature.  The  central,  longitudinal  channel  axis  at  the  level  of  the  bed  has 
constant  mean  slope  Sc,  and  describes  a  helix  in  space  which  traces  a  circle 
of  radius  rc  when  projected  onto  a  horizontal  plane.  The  flow  is  conveniently 
described  in  cylindrical  coordinates:  the  vertical  z  axis  passes  through  the 
curvature  center  of  the  channel  and  is  positive  in  the  direction  opposite  to 
gravity;  in  planes  perpendicular  to  the  z  axis,  locations  are  specified  by 
radial  distance  from  the  z  axis,  r,  and  polar  angle,  e,  as  shown  in  figure  1. 
In  order  for  the  radial  slopes  of  the  bed  and  water  surface  to  be  constant 
along  the  channel,  the  local  streamwise  slope,  S(r),  of  both  must  be 
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s  =  s. 


The  flow  is  treated  as  uniform  in  the  sense  that  its  properties  are 
invariant  along  any  helix  with  constant  radius  r  and  slope  S.  The  analysis 
will  be  restricted  further  to  a  central  region  of  the  channel,  delineated  in 
figure  1,  throughout  which  the  vertical  velocity  is  much  smaller  than  the 
characteristic  velocities  in  the  r  and  e  directions.  The  channel  slope  Sc 
will  be  limited  to  small  values,  so  that  forces  and  velocities  parallel  to  the 
underlying  bed-surface  helices  may  be  taken  to  be  equal  to  those  along  the  0 
coordinate.  Finally,  the  restriction  d/r  <<  1  will  be  imposed,  for  reasons 
that  become  apparent  in  the  next  section. 

Vertical  distribution  of  radial  shear  stress.  Calculation  of  the  radial 
shear  stress  at  any  elevation  requires,  first,  that  the  radial  water-surface 
slope  and  associated  pressure  gradient  be  determined,  for  use  in  calculation 
of  the  vertical  distribution  of  radial  shear  stress  from  the  radial  momentum 
equation.  The  radial  water-surface  slope  will  be  calculated  from  a 
simplified,  by  means  of  an  ordering  analysis,  formulation  of  the  conservation 
of  flux  of  moment-of-momentum  for  the  control  volume  shown  in  figure  1,  which 
extends  over  the  whole  flow  depth  and  has  base  dimensions  &r  and  as  =  r  A6. 
For  this  control  volume  equation  of  moment-of-momentum  about  an  axis  r  = 
constant  at  the  bed  surface  is 


1  „  12  1  ,  1 

3p  V  1  3  r  j2  r  .  1 
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pressure; 


where,  in  addition  to  the  quantities  defined  in  figure  1,  p  = 

n  =  ~p;  p  =  fluid  density;  u(r,n),  v(r,n),  and  w(r,n)  =  velocities  in  r, 

o,  and  i  directions,  respectively;  and  Tp2(r,n)  =  T2r(r,n)  =  shear  stress 

acting  on  surfaces  normal  to  r  and  z  axes,  respectively.  The  fourth  term  in 

(2)  is  zero  for  uniform  flow.  The  remaining  terms  will  be  ordered  by  taking  v 

-  0(V),  where  V(r)  =  depth-averaged  flow  velocity;  u  =  0(u(r,l)  =  Um);  and 

the  z-direction  velocity  w  =  0(w  =  W  ).  Relative  to  the  second  term,  the 

J  '  max  nr 

Um  r  Um  Wm 

fifth  and  sixth  terms  are  of  order  0(-7j— )  and  °(q  y~)»  respecti vely.  Yen 

(1965)  concluded  from  his  measurements  of  flow  in  curved  channels  that 

U 

•y-  =  O(-p).  This  result  is  suggested  also  by  the  analysis  of  Zimmermann  and 
Kennedy  (1978,  Eq.  7),  which  shows  the  ratio  of  radial  to  streamwise 

components  of  bed  shear  stress  to  be  O(-p).  If  z  =  0(d),  the  continuity 

equation. 


u 
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(3) 


in  which  is  zero  for  the  uniform  flow  being  considered,  requires 
W  3  S  U  W  2 

that  jj-  =  0 (“) .  Because  ^  =  0(~),  it  follows  that  =  O(^)  .  It  is 

concluded  then  that  the  fifth  and  sixth  terms  of  (2)  are  both  0(-^)  relative  to 

the  centrifugal-force  (second)  term,  and  can  be  dropped. 

The  shear-stress  (third)  term  of  (2)  may  be  ordered  by  utilizing  the 

equality  of  shear  stresses  and  the  Boussinesq  shear-stress  relation  for 

turbulent  flow,  and  treating  the  eddy  viscosity,  eq,  as  constant  over  the 

depth: 
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1  Trz  dn  =  l  Tzr  dri  =  p  /  l  §7  dr>  =  p  T-3 

The  eddy  viscosity  may  be  expressed  as  a  product  of  the  shear  velocity,  u*, 
and  depth;  therefore, 

1 

f  t  dn  =  0(apu*U)  (5) 

0 

where  a  =  eQ/u* d  (a  =  0.079,  according  to  Hinze  (1975)).  From  (5)  it  follows 
that 


in  which  f  =  Darcy-Weisbach  friction  factor.  Both  a  and  /f/8  are  0(10  1),  and 

therefore  the  third  term  of  (2)  is  two  orders  of  magnitude  smaller  than  the 

2 

second  and  may  be  disregarded.  Because  p  =  0(pV  ),  the  first  and  second 
terms  are  of  the  same  order.  Incorporation  of  the  simplifications  resulting 
from  the  foregoing  ordering  analysis  reduces  (2)  to 


/  i  dn 
0  8r 
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In  the  central  region,  where  ~  «  1  and  «  1,  it  is  reasonable  to 


assume  that  the  vertical  distribution  of  p  is  hydrostatic: 
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where  g  =  gravitational  acceleration  and  H  is  defined  in  figure  1.  It  has 
been  demonstrated  by  Yen  (1965)  that  the  deviation  of  the  pressure  from  the 
hydrostatic  distribution  is  O(-p)  or  smaller  in  even  moderately  curved  open- 
channel  flow.  The  primary-flow  velocity,  v(r,n),  will  be  expressed  by  the 
power  law. 


v  _  n+1  1/n 


where  1/n  =  exponent,  which  is  related  to  the  Darcy-Weisbach  friction  factor 


-  =  -  /f78 
n  k 


where  k  -  Karman's  constant.  The  background  of  this  relation  is  reviewed  by 
Zimmermann  and  Kennedy  (1978).  Karim  (1981)  examined  (10)  critically,  veri¬ 
fied  it  with  laboratory  data,  and  formulated  the  dependence  of  <  on  sediment 
concentration;  this  refinement  will  not  be  included  in  the  present  analysis. 
Substitution  of  (8)  and  (9)  into  (7)  yields 


H'  (r)  -  Hill- 
'  '  n  rg 


By  means  of  an  ordering  analysis  similar  to  the  one  developed  above  and 
guided  in  some  measure  by  his  experimental  results.  Yen  (1965)  simplified  the 
radial  momentum  equation  for  curved  open-channel  flow  to 
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It  is  noteworthy  that  xzr  makes  a  first-order  contribution  to  the  radial- 
momentum  relation,  but  the  corresponding  vertical  shear  stress,  x  ,  makes 
only  a  higher-order  contribution  to  the  moment-of-momentum  equation  if  the 
moment  axis  is  taken  at  bed  level  to  avoid  inclusion  of  the  bed  shear 
stress.  This  is,  in  fact,  the  motivation  for  utilizing  the  moment  relation: 
it  avoids  specification  of  x  at  the  bed  in  the  determination  of  H'. 
Substitution  of  (9)  into  (12),  integration  of  the  resulting  expression  from 


arbitrary 


n  to  0  =  1,  and  application  of  the  boundary 


condition  x  (1)  =  0  yields 


2  2 

.(r,n)  =  pgd  [H'(n-l)  -  ~  ^+2 )  “  1 ) J 


Traverse  bed  profile  and  depth  distribution.  Equilibrium  of  the  trans¬ 
verse  bed  profile,  h(r,e)  in  figure  1,  and  of  the  depth,  d(r),  are  attained 
when  the  radial-plane  forces  acting  on  the  moving,  bed-load  particles  sum  to 
zero.  Bed-load  movement  will  be  treated  as  occurring  in  a  layer  of  thickness 
y b ,  as  shown  in  figure  1.  The  shear  forces  exerted  on  this  agitated,  somewhat 
dilated,  moving  layer  are  in  reality  diffuse,  "seepage"  forces  caused  by  the 
flow  within  the  layer's  intersticies,  and  any  net  force,  however  small,  in  the 
radial  direction  will  produce  transverse  motion  of  the  bed-load  particles. 
Therefore,  radial  equilibrium  will  be  reached  when  the  local  bed  inclination, 
B(r),  is  such  that 


Tor  =  Tzr(°)  =  yb(1“P)  ap  9  sin  B 


(14) 
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where  p  =  porosity  of  the  bed-layer;  A p  =  p$-p;  and  =  density  of  bed  parti¬ 
cles.  In  the  development  of  his  detailed,  computer- simulation  model  of  sedi¬ 
ment  transport  in  streams,  Karim  (1981)  concluded  from  inferential  evidence 
that 


(15) 


where  D50  =  median  bed-material  size,  u*(r)  =  local  shear 

velocity  =  V/f/8;  and  u*c  =  critical  shear  velocity  for  incipient  particle 
motion.  In  terms  of  the  Shields  parameter,  0,  u*c  may  be  expressed 


(16) 


The  traverse  bed  profile  may  be  calculated  by  neglecting  the  effect  of  H' 


on  d(r);  then 


•  „  dd 

sin  0  =  dr 


The  velocity  V  to  be  used  in  (18)  is  obtained  by  incorporating  (1)  into  the 
local  expression  of  the  Darcy-Wei sbach  friction  factor. 


V (r)  =  /8 


/ss 


where  tq6  =  longitudinal  shear  stress  acting  on  the  bed;  and  6(r)  =  bed-shear- 
stress  reduction  factor  defined  by 


Tnfl(r)  =  <SSPgd(r) 


which  takes  into  account  the  transport  of  primary-flow  momentum  out  of  the 
central  region  to  the  vicinity  of  the  outer  bank,  where  it  is  balanced  by  bank 
shear.  An  analysis  of  6  is  developed  below.  Substitution  of  (19),  (20),  and 
(1)  into  (18)  and  integration  of  the  resulting  expression  for  ^j~-  yields 


1  /8e  1  +  /f  AW 


1  1  r  1  1  -  / 


/d  /rc  /?  /rc  (1_p)  i  +  2/7  '  f  g  d50 


where  the  subscript  c  denotes  the  centerline  values  used  in  setting  the  inte¬ 
gration  constant.  Elimination  of  6SC  from  (22)  by  means  of  (20)  and  replacing 
V  by  its  section  averaged  value  for  the  whole  flow,  7, to  facilitate 


verification,  leads  to 


(23) 
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M-  =  1  -  (1  -  /~f-)  !£P^t  ~  Fn 


r.'  (1-p) 


1  +  2/f 


where  Fn  =  V/yg^-D 


D  ~  ▼  a  p  “50 

Utilization  of  the  mean  velocity  for  the  whole  cross  section  in  the 


calculation  of  d('*)  from  (23),  or  in  calculating  an  average  value  of  Sy  from 
(18),  by  replacing  FD  with  neglects  the  effects  of  the  variation  of  V 
across  the  channel,  but  nevertheless  yields  satisfactory  results,  as  is 
demonstrated  in  Section  III. 

Equations  20  and  23  give  the  radial  distribution  of  mean  velocity  for 
uniform  flow.  In  practical  applications,  it  generally  suffices  to  take 
6  =  1. 

Vertical  distribution  of  transverse  velocity.  Calculation  of  the  radial- 
plane  velocity  will  incorporate  the  following  assumptions: 


1.  The  primary-flow  shear  stress,  t  ,  is  linearly  distributed 
9t  29 

and  makes  a  negligible  contribution  to  the  streamwise  force  balance: 


Tz0(r,n)  =  T  oe^1 2"11^  =  <5P9dsO-ri) 


(24) 


2.  The  eddy  viscosity  is  isotropic  and  is  given  by 


e(r»n)  =  7  26 


p  9_v 
9n 


(25) 


3.  Because  of  the  isotropy  of  e,  the  radial  velocity  and  shear  stress 
are  related  by 
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Substitution  of  (9)  and  (24)  into  (25)  yields 

„  ,2  _  _  n-1 

<5gd  Sr  2  — - 
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e - rV  n+T 


which,  when  substituted  along  with  (13)  into  (26),  leads  to 

l  —  l~n 

il  _  I £_  nil  ’r  r  h'  n  —  fo+1 ^  nn  -n  n  ^  -i 
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For  steady,  uniform  flow,  u  must  satisfy 


/  u(n)  dn  =  0 
0 


There  is  no  assurance  that  (28)  will  satisfy  this  requirement  if  H *  given  by 
(11)  is  utilized,  because  of  errors  inherent  in  the  Boussinesq  eddy-viscosity 
model  and  other  assumptions  that  have  been  made  in  the  derivation  of  the 
relation  for  u.  Therefore,  the  integral  of  (29)  will  be  evaluated  for  arbi¬ 
trary  H',  denoted  by  H’u  and  expressed  as 


H;<r>5  T^7g 


where  T  =  H^/H '  and  H'  is  given  by  (11),  Substituion  of  (28)  and  (30)  into 
(29),  utilization  of  the  expansion 
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(31) 


and  term-by-term  integration  of  the  resulting  series  yields 
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Equation  29  is  satisfied  if  T  is  given  by 


T(n)  =  -  |ntl)  E  [  3 
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]  (33) 


Incorporation  of  (17),  (20),  and  (33)  into  (32)  gives 
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Equation  34  is  portrayed  in  figure  2  for  four  values  of  n  which  span  the  range 
from  very  rough  (n  =  2.5;  f  =  0.16)  to  relatively  smooth  (n  =  10,  f  =  0.01) 
channels.  It  is  noteworthy  that  for  all  but  very  low  values  of  n,  the 
velocity  profiles  are  nearly  linear  except  near  the  bed,  with  u  =  0  at  about 
mid-depth. 

A  remark  on  H1  and  H'u.  In  the  foregoing  analysis,  two  expressions  were 
derived  for  the  transverse  slope  of  the  water  surface:  (11),  and  (30)  and 
(33).  Corresponding  to  each  of  these  is  a  different  value  of  Tzp(0),  the 
radial  component  of  the  bed  shear  stress.  Their  ratio  T  =  H ' u /H '  given  by 


A14 


(33)  has  a  nearly  constant  value  of  0.9  for  2  <  n  <  8.  In  view  of  the  near 


equality  of  H'u  and  H1,  why  was  it  necessary  to  utilize  different  values  in 


the  formulations  of  tzr  and  of  u(n)  ?  The  problem  is  one  of  sensitivity,  as 


will  now  be  demonstrated.  Equation  13  gives 


«(r)  =  tzr(r,0)  =  pgd  [75^%-  H‘] 


which  together  with  (11)  for  H'  shows  that  t  is  the  difference  between  two 

or 


small  quantities  multiplied  by  a  large  one,  (pgd).  Therefore,  small  errors  in 


the  expression  for  the  transverse  water  slope  produce  large  errors  in  x  .  For 


example,  the  ratio  of  t  given  by  (35)  with  H'  replaced  by  H'u  obtained  from 


(30)  and  (33),  to  xQr  yielded  by  (35)  and  (11)  varies  widely  with  n,  from 


about  0.6  for  n  =  2  to  nearly  0.2  for  n  =  8.  Because  the  radial  bed  slope  and 


thus  also  the  bed  profile  depend  directly  on  xQp  ,  as  is  shown  by  (14),  it  is 


important  in  their  derivation  to  have  an  accurate  estimate  of  x  --one  whose 

or 


calculation  avoids  use  of  such  artifices  as  the  Boussinesq  relation,  and 


instead  directly  utilizes  a  mechanics  principle  such  as  conservation  of  moment 


of  momentum.  The  effect  of  the  radial  water-surface  slope  on  u(n)  calculated 


from  the  Boussinesq  relation  may  be  examined  by  substituting  (26)  into  (12) 


and  treating  e  as  constant,  say  eQ  ,  which  results  in 


32u  /V2  m,  v 

{v'*) 


Equation  36  shows  that  L_u  ,  and  hence  u(n)  »  also  Is  very  sensitive  to  the 

3nZ 

small  difference  between  two  nearly  equal  quantities  multiplied  by  a  large 


one.  Accordingly,  only  a  very  small  adjustment  in  the  radial  water-surface 


/.  .  \ 
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slope  is  required  to  compensate  for  the  effects  on  u  of  the  assumptions  made 
in  its  calculation,  and  thereby  to  permit  u  to  satisfy  the  continuity 
requirement,  (29).  But  this  small  adjustment  in  H'--amounting  to  only  about 
10  percent,  as  just  noted--has  a  major  effect  on  TQr  ,  as  (35)  demonstrates. 

Secondary-flow  effects  on  t  and  v.  In  a  laterally  nonuniform  curved 

flow,  the  secondary  current  produces  a  net  radial  transport  of  streamwise 

(e-direction)  momentum  out  of  the  central  region,  which  in  turn 

reduces  tq  and  in  so  doing  modifies  v(r,n)(or  n).  It  was  anticipation  of  this 

effect  which  prompted  incorporation  of  the  factor  6  (r)  into  the  shear-stress 

expressions,  (20)  and  (21).  Calculation  of  6  and  the  secondary- flow  effect 

on  v(n)  proceeds  from  the  e  -direction  momentum  equation  with  j  neglected, 

10 


„  3v  3_v  ,  3  v  jjv 
u3r  v3s  w3z  r 
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Multiplication  of  (3)  by  v,  addition  of  the  result  to  (37),  integration  of  the 
new  relation  from  n  to  1,  and  imposition  of  the  boundary  condition 
TzQ(r,l,0)  =  0  gives 


^  r  ,  3 (uv) 
Tze  "  pd  {-  ; 


dn  -  —  /  uv  dn 
n 


1  ,  2 
,  3  v 

1  FT  * 
n 


d  ^wv^n=l  '  wv^  +  9  S(l-n)}. 


where  S  =  -  r-r.  The  z  velocity,  w,  is  evaluated  from  the  continuity  equation, 

d  S 

(3),  the  relation  for  u,  (32),  and  the  identities, 


and,  similarly. 


9n  _  S 
as  cf* 


The  result  is 


/  \  _  rd  ,  dd  . /3  dV  1  d6.\-i  r  j 

w(r,n)  “  “  [7  +  df  +  d{-  d?  "  6  d?)]  ^  u  ^ 

r3H  dd  / «  \  1  c 

+  [?F  ‘  d?  (1_n)]  u  "  Sv* 


Substitution  of  (9),  (34),  and  (41)  for  v,  u,  and  w  along  with 


2  2 
av  s  av 


as  d  an 


into  (38)  yields 


•  (l-n )  -  646  f  n(n+l){[^*  2  7  *  4  y  37  -  £  fj  '  G(n,„)  nn  dn 


,  rdd  ,  d  .  -  d  dV  d  d6n  n  .  r,  ,  ,  , 
*  [d?  +  7  +  3  V  dF  •  s  d?]  n  '  G(n’'')  dn> 


where  G  is  defined  by  (34). 

The  first  terms  on  the  right-hand  side  of  (43)  is  the  linear  shear-stress 
distribution,  and  the  second  term  expresses  the  shear-stress  reduction  due  to 
the  transverse  gradient  of  lateral  flux  of  streamwise  momentum.  In  the  calcu- 

<\  jp 

lation  of  t  .  it  is  assumed  that  the  f—  term  is  negligible  in  comparison  to 
Z0  3  r 

the  other  derivative  terms  in  brackets.  Substitution  of  (18)  for  4^  and  of 


(20)  for  V  then  permits  calculation  of  the  shear  stress  at  any  location.  The 
bed-shear-stress  reduction  factor,  6,  introduced  in  (20)  and  (21),  is  obtained 
by  letting  n  =  0  in  (43).  To  gain  some  idea  of  its  magnitude,  a  representa¬ 
tive  constant  value  of  6  for  the  whole  central  region,  say  6,  will  be  deter¬ 
mined.  For  this  calculation  it  is  appropriate  to  replace  V  and  d  by  their 
section-averaged  values,  V  and  d,  after  carrying  out  the  substitutions  and 
taking  the  derivatives  in  (43);  and  to  take  r  =  rc.  The  result  is 

t  -  1 

6  =  =[1  +  192  7  n  (n+1 )  ST  /  G(n,n)n1/n  dr,]"1  (44) 

pgba  rc  i  q 

where  Sy  is  obtained  from  (18)  by  replacing  d  and  V  by  d  and  V.  The  integral 
of  (44)  was  evaluated  numerically,  with  the  result  shown  in  figure  3.  Values 
of  1  for  some  field  and  laboratory  flows  are  presented  in  the  next  section. 

The  effect  of  the  secondary  flow  on  the  primary-flow  velocity 
distribution  may  be  estimated  by  substituting  (17)  into  (20)  and  replacing 
<5  by  J.  If  V  is  considered  to  be  constant,  n  is  increased,  as  1//6  ,  which 
corresponds  to  the  velocity  becoming  blunter.  This  is  the  observed  effect  of 
secondary  flow  or  v(n)  (Falcon  1978). 

III.  VERIFICATION 

Data  utilized  in  the  verifications  reported  here  are  summarized  in  table 
1.  Falcon  (1978)  presents  additional  comparisons  of  measured  and  computed 
quantities. 

Bed  topography.  Zimmermann  (1974)  and  Zimmermann  and  Kennedy  (1978) 
reported  the  results  of  experiments  conducted  in  three,  concentric,  60-cm 
wide,  circular-plan-form  flumes  with  a  central  angle  that  approached  2ir .  Two 


different  sediments,  with  median  diameters  of  0.21  mm  and  0.55  mm,  were 
used.  Longitudinal  bed  profiles  were  measured  at  11  different  radii,  and  the 
transverse  bed  profiles  obtained  from  them  for  numerous  cross-sections  in  the 
reaches  of  fully  developed  flow  were  plotted  and  averaged.  The  transverse 
profiles  were  found  to  be  slightly  convex  upward,  as  illustrated  in  figure 
4.  Mean  transverse  bed  slopes,  averaged  across  numerous  sections  for  each 
run,  were  then  computed.  Figure  5  shows  the  transverse  slopes,  Sy  ,  so  deter¬ 
mined  plotted  in  the  format  of  (18)  based  on  cross-section-averaged 
properties.  Excellent  agreement  between  measured  and  computed  values  is 
obtained  if  ^  ^  =  1.3.  If  the  limiting  value  (for  fully  turbulent  boundary 
layers)  of  the  Shields  parameter,  e  =  0.06,  is  adopted,  the  resulting  porosity 
is  p  =  0.47,  a  not  unreasonable  value  for  the  agitated,  dilated,  moving  bed¬ 
load  particles.  The  computed  profile  for  Zimmermann's  (1974)  Run  No.  RII-13 
shown  in  figure  4  was  obtained  from  (23),  using  these  values  of  e  and  p.  The 
centerline  depth,  dc  =  9.66  cm  utilized  in  computing  the  profile  was  obtained 
by  equating  the  reported  mean  depth,  10.1  cm,  to  the  mean  depth  d  calculated 
by  integration  of  d  given  by  (23)  across  the  channel  width: 


f-  =  1  -  2*  +  2*  2  +  | 


2„  3/2  3/2, 

’♦  Hr0  -  r,  ) 


/r~  (r  -  r . ) 
c  '  o  i ' 


(45) 


where  $  =  Fn  TT9,n'i  ~~~~i  ,  and  r^  and  rQ  =  radii  of  the  inner  and  outer 


banks,  respectively.  Calculation  of  the  relation  between  d  and  dc  in  this  way 


is  consistent  with  the  measurement  procedure  that  was  used.  Falcon  (1978) 
describes  calculation  of  "d  in  a  way  that  is  consistent  with  conservation  of 
bed-material  volume  in  a  curved  channel.  The  friction  factor  utilized,  f  = 
0.165  for  this  flow,  is  the  value  for  the  bed  section  obtained  from  the  side- 
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wal 1 -correct  ion  procedure  (Vanoni  1976,  p.  152).  The  measured  and  computed 
profiles  are  in  excellent  agreement.  The  bed-shear-stress  reduction  factor 
obtained  from  (44)  for  this  flow  is  6  =  0.43,  which  shows  that  in  this 
relatively  narrow  channel  the  secondary  current  produced  a  major  reduction  in 
the  bed  shear. 

Figure  6  compares  measured  and  computed  transverse  bed  profiles  for  a 
Missouri  River  section  (Falcon  1978).  In  computing  the  profile,  d  was  taken 
to  be  equal  to  dc,  because  of  the  difficulty  of  determining  rc  for  a  wide 
natural  stream,  and  of  the  insensitivity  of  dc  to  rc  (see  (45)).  Included  in 
figure  6  is  the  mean  bed  profile  given  by  (18),  which  also  is  seen  to  give 
quite  good  results.  Equation  44  gives  6  =  0.97  for  this  relatively  wide, 
shallow  flow,  demonstrating  the  bed  shear  stress  at  any  r  in  the  central 
region  of  this  natural  channel  is  nearly  equal  the  local  value  of  pgdS.  Other 
comparisons  of  measured  and  computed  bed  profiles  yielded  conformities  as  good 
as  those  demonstrated  in  figures  4  and  6.  In  evaluating  data  from  natural 
streams,  in  which  the  flow  is  seldom  steady  for  appreciable  periods,  there  is 
always  uncertainty  about  the  equilibrium  of  the  bed  topography.  Moreover,  the 
bed-material  size  often  varies  widely  across  a  section,  often  by  a  factor  of 
five  or  more.  In  view  of  these  difficulties,  it  is  suggested  that  in  the 
calculation  of  bed  topography  by  means  of  (18)  and  (23),  averaged  (across 
channel  sections,  and  along  subreaches  that  are  sufficiently  short  that  rc  is 
practically  constant)  values  of  d  and  V  be  used,  and  that  the  median  diameter 
of  the  material  that  can  be  moved  by  the  flow  be  utilized  for  D^q. 
Furthermore,  for  most  natural-stream  situations,  the  refinement  given  by  (23) 
is  probably  not  justified;  a  straight-line  profile  with  slope  Sy  given  by  flHl 
and  passing  through  d  =  d  at  r  =  rc  is  generally  as  accurate  as  the  *  ^ 
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data  warrant,  and  perhaps  within  the  reproducibility  of  bed  topography  of 
natural  streams  with  their  vagaries  of  discharge  and  bed-sediinent 
characteristics. 

Velocity  distributions.  It  is  very  difficult  to  obtain  reliable  data  on 
u(o,r)  in  erodible-bed  channels,  because  of  the  small  values  of  the 
secondary-flow  velocities,  and  of  the  problems  posed  by  the  moving  sediment 
and  the  continuous  bed  changes  attendant  to  migration  of  bed  forms. 
Therefore,  the  two  measured  profiles  obtained  by  Kikkawa,  Ikeda,  and  Kitagawa 
(1976)  from  uniform  flow  in  a  circular-plan-form,  rigid  channel  were  utilized 
in  the  verification  of  (34),  with  the  results  shown  in  figure  7.  Kikkawa  et 
al  (1976)  developed  an  analytic  model  for  u(n)  ,  which  can  be  seen  in  their 
paper  also  to  yield  generally  satisfactory  results  except  near  the  bed,  where 
it  does  not  satisfy  the  no-slip  condition.  Comparisons  presented  by  Falcon 
(1979)  of  (34)  with  the  rigid,  sinuous-channel  data  on  transverse-velocities 
reported  by  Yen  (1965)  also  demonstrate  very  satisfactory  agreement. 

The  transverse  distributions  of  V,  the  depth-averaged  streamwise 
velocity,  in  erodible-bed  channels  are  somewhat  easier  to  measure  than  the 
radial-velocity  distributions.  Velocity  data  obtained  by  Onishi  (1972)  at  the 
apex  cross  sections  in  two  of  his  rigid-bank,  erodible-bed,  meandering-channel 
flows  were  used  to  validate  the  distribution  of  V  obtained  from  (20)  and 
(23).  The  average  friction  factors,  7,  used  in  the  computations  were  obtained 
from  the  reported  mean  values  of  velocity,  depth  and  slope  for  the  flows,  and 
the  Darcy-Weisbach  relation  in  the  form 

7  =  MS  (46) 

r 
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Note  that  computation  of  7  from  (46),  which  assumes  6=1,  and  of  the 
corresponding  n  from  (17),  computation  of  J  from  (44)  for  this  n,  and  use  of 
fi"  and  f  in  (20)  to  determine  V(r)  is  slightly  inconsistent.  Falcon  recommends 
iteration  between  (20)  with  V  =  V  ,  (44),  and  (17)  to  obtain  consistent  values 
of  n,  f,  and  T  .  However,  the  convergence  is  quite  rapid  and  the  effect  on 
the  computed  V ( r )  or  u(n)  is  not  great.  The  computed  and  measured  distribu¬ 
tions  of  V(r)  shown  in  figure  8  agree  quite  well  except  near  the  banks,  where 
V  must  tend  to  zero. 

IV.  CONCLUDING  REMARKS 

It  must  be  borne  in  mind  that  the  model  developed  here  is  strictly  valid 
only  for  uniform,  curved-channel  flows.  However,  the  available  experimental 
data  on  flows  in  strongly  curved  channels  (Zimmermann  1974,  Zimmermann  and 
Kennedy  1978,  Odgaard  and  Kennedy  1982)  indicate  that  they  are  characterized 
by  relatively  small  phase  shifts  or  lag  distances  between  local  secondary-flow 
properties  or  bed  topography  and  local  channel  curvature.  Therefore, 

application  of  the  present  model  utilizing  local  channel  and  flow 
characteristics  in  nonuniform  flows,  as  was  done  in  the  foregoing  comparison 

with  Onishi's  (1972)  data,  will  generally  yield  satisfactory  results. 

However,  in  the  case  of  flow  in  weakly  meandering  sinuous  channels,  as 
investigated  by  Gottlieb  (1976)  and  Falcon  (1979),  the  phase  shift  between 

local  channel  curvature  and  secondary-flow  strength  approaches  tt/2. 

The  analytical  model  developed  here  is  valid  only  for  the  central  regions 
of  curved-channel  flows,  which  generally  extend  to  about  one  local  depth  from 
the  bank.  In  the  near-bank  regions,  the  flow  becomes  strongly  three- 
dimensional  and  heavily  influenced  by  local  bank  characteristics  (erodibility. 


slope,  roughness,  etc).  Analysis  of  the  flow  and  sediment  transport  in  these 
regions  is  correspondingly  more  difficult  than  for  the  central  region,  and 
likely  must  await  availability  of  better  experimental  data  for  its  guidance. 
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FIGURE  CAPTIONS 

1.  Definition  sketch  for  flow  in  an  alluvial-channel  bend. 

2.  Distribution  of  radial-plane,  secondary-flow  velocity  given  by  (34). 

1  £ 

3.  Result  of  the  numerical  evaluation  of  f  G(n,n)n  dn  ,  where  G  is 

given  by  (34).  0 

4.  Transverse  bed  profiles  measured  in  Zimmermann's  (1974)  Run  RII-13, 
and  computed  from  (23). 

5.  Comparison  of  (18)  and  Zimmermann  and  Kennedy's  (1978)  measured 
transverse  bed  slopes. 

6.  Transverse  bed  profile  measured  in  Missouri  River  (Falcon  1978)  and, 

those  computed  from  (18)  ( - )  and  (23)  ( - )  . 

7.  Secondary- flow  velocity  profiles  measured  by  Kikkawa  et  al  (1976)  and 
computed  from  (34). 

8.  Comparison  of  Onishi's  (1972)  measured  transverse  distributions  of 
depth-averaged  velocity  and  those  computed  from  (20)  and  (23). 


Figure  2.  Distribution  of  radial-plane,  secondary-flow  velocity  given  by  (34). 


0  1  + 

COMPARISON  or  MEASURED  AND  COMPUTED  TSANS  VLKSE  BED  SLOPES  IN  ALLUVIAL  CHANNELS 

(Data:  Zimmerman*  and  Kennedy  1978) 

Figure  5.  Comparfson  of  (18)  and  Zfmroermann  and  Kennedy's  (1978) 
measured  transverse  bed  slopes. 


igure  6.  Transverse  bed  profile  measured  in  Missouri  River  (Falcon 
those  computed  from  (18)  ( - }  and  (23)  ( - ). 
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APPENDIX  B:  LISTING  OF  COMPUTER  PROGRAM  PR-SEG6 


AND  INPUT-OUTPUT  SAMPLES 

MAIN  PROGRAM  PR-SEG6 


*“**»**»♦  MAIN  PROGRAM  PK-SEG6  **.***.*******.****** 

THIS  1*5  A  GENERAL  PROGRAM  FOR  A  RIVER  REACH  COMPOSED  * 
‘"’F  SEVERAL  SEGMENTED  BENDS  WITH  POSITIVE  OR  NEGATIVE  * 
RADII  OF  CURVATURE.  THIS  PROGRAM  SIMULTANEOUSLY  * 
SOLVES  THE  CONTINUITY  AND  MOMENTUM  EQUATIONS  BY  * 
ITERATING  BETWEEN  THEM.  EACH  NEW  SEGMENT  INPUT  DATA  * 
MUST  BE  READ  AT  THE  FIRST  SECTION  OF  EACH  SF  G  ME  NT •  * 


C  «  * 


*********** 


**************** 


C  A  lq 


subroutines:  * 

FALL  :  DETFRMINES  THE  F -P AR AME T ERS *  WHICH  DEPEND  ON  THE* 
FOLLOWING:  * 

FI  =  <  ST,  RC,  R,  H  )  * 

F 2  -  (  ST,  RC*  R,  Gl,  G2,  G3,  H  )  * 

F?  =  <  ST,  RC*  V,  R,  Gl,  G2 ,  G3»  H,  F,  VBA  R  )  * 

F  4  =  <  ST,  R,  H  )  * 

GALL  :  DETERMINES  THE  G-P AR AME T E R S ,  WHICH  DEPEND  ON:  * 

G  =  C  RHOS,  BETA,  ALPHA  )  * 

Gl  =  <  F*  BETA  )  * 

G2  =  (  F  )  * 

G 3  =  (  VBAR,  RHOS,  DEO,  F,  BETA,  ALPHA,  THETAC*  * 
POR ,  G  )  * 

INTGRL  :  EVALUATES  THE  INTEGRAL:  * 

SUM  =  I  NT (  R  *  S  OR  T  (  C A  *  R  ♦  B )  /  ( R  ♦  D)  )  )  * 

FROM  R(J-l)  TO  R ( J )  * 

C  A  LQ  :  DETFRMINES  THE  DISCHARGE  OVER  THE  CROSS  * 

SECTION  given:  * 

C  M,  V,  HL*  R  -  GETS  VA,  VQ,  AT,  QT  )  * 

***************************** 
variables:  * 

l» H I  =  TRANSVERSE  SHIFT  VFLOCITY  (FT/SEC)  * 

UBIM1  =  VALUES  FOR  UB AR  AT  THE  PREVIOUS  SECTION.  USED  * 
WHEN  CHANGING  ORIENTATION  OR  WHEN  USED  AS  AN  * 

APPROXIMATION  TO  THE  CURRENT  SECTION'S  VALUES  * 

<KOPTl=l>  * 

TUP  =  TEMFORARY  STORAGF  op  U3AR  WHEN  CHANGING  * 

ORIENTATION  * 

UB M I  =  UPDATED  TRANSVERSE  SHIFT  VELOCITY  (FT/SEC)  * 

VI  =  STRCAHUISE  VELOCITY  (FT/SEC)  * 

V I M 1  =  PREVIOUS  VALUE  OF  VI  AT  (  I  -  1  )  * 

TE  M V  =  TEMPORARY  STORAGE  OF  V  WHEN  CHANGING  ORIENTATION* 
OR  VALUES  OF  ETA-MODIFIED  V  * 

V  M I  =  UPDATED  STREAMWISE  VELOCITY  (FT/SEC)  * 

UI  =  LOCAL  SECONDARY  FLOW  VE LOC I T Y  —  MOD  I F I E C  (FT/SEC)* 
HL  I V  =  LOCAL  PLOW  DEPTH  AT  SECTION  I  (FT)  * 

NONPIMFNSIONALIZED  BY  THE  VERTICAL  LENGTH - H  * 

HL I H  =  LOCAL  cLOw  DEPTH  AT  SECTION  I  <PT>  * 

NONE.  IMENSIONALIZED  BY  THE  HORIZONTAL  LENGTH  — -W  * 
LIM1V  =  LOCAL  PLOW  DEPTH  AT  PREVIOUS  SECTION  (1-1)  FT  * 


URN  I  = 
VI  - 
VI  Ml  = 
TE  M  V  = 

VN  I  = 
UI  = 
HL  I  V  = 

HL  I  H  = 


HL  I M  IV  =  LOCAL  PLOW  DEPTH  AT  PREVIOUS  SECTION  ( 
HLIM1H  =  "  "  "  *  »  " 

THLV  =  TEMPORARY  STORAGE  OF  HL V  WHEN  CHANGING 
ORIENTATION 

THLH  =  «  "  "  HLH  "  " 
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n  o 


C 

C 

C 

c 


c 

c 

c 

c 

c 

r 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c  * 

c 

c 

c 

c 

c 

c 

c 

c 

c 

r 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

r 

c 


SI  = 
SIM1  = 

TEMPS  = 

uc  = 
uo  = 

STI  = 
R  = 
TEMPR  = 
VA  = 
VQ  = 
V  = 
VS  = 
U  - 
HL  “ 
ILOC  = 

1ST  = 


ISEG  = 
NIST  = 
ISTD  = 

IKE  E  P  = 


H  = 
HV  = 
HH  = 
PC  = 


R  n  r 

$1*  S  2  = 

SC  = 
R  C  D  = 
OS  r 
F  = 


VBAR  - 

R  C  L  - 

U  = 


LOCAL  STPEAMWISE  POSITION  AT  SECTION  I  (FT  )  * 
LOCAL  STREAMWISE  POSITION  ST  PREVIOUS  SECTION 

(I  -  1)  <  F  T  ) 

TEMPORARY  S  STORAGE  FOR  REVERSING  ORIENTATION 


CENTERLINE  SECONDARY  FLOW  VELOCITY  (FT/SEC)  * 
INITIAL  CONDITION  FOR  THE  SECONDARY  FLOW  * 
INTEGRATION  * 
CUMULATIVE  TRANSVERSE  BED  SLOPE  * 
RADIAL  COORDINATE  FROM  THE  CENTERLINE  (FT)  * 
TEMPORARY  R  STORAGE  FOR  REVERSING  ORIENTATION  * 
LOCAL  FLOW  AREA  <SQ  FT)  * 
LOCAL  PLOW  DISCHARGE  CCFS)  * 
STORAGE  ARRAY  FOR  OUTPUT  * 

n  n  n  n  , 
H  n  n  n  * 
n  n  n  n  4 


ARRAY  DENOTING  THE  SECTION  NUMBER  WHERE  EACH  NEW* 


SEGMENT  BEGINS  * 

ARRAY  DENOTING  THE  SECTION  NUMBER  WHERE  THT  * 

VELOCITIES  AND  DEPTHS  ARE  KEPT  UNTIL  FINAL  * 

TABULATION  MAXIMUM  OF  6  VALUES  STOR ED -T NLE T  ♦  * 

A  OTHERS*  OUTLET  * 

VARIABLE  FOP  THE  ILOC  ARRAY  * 

VARIABLE  FOR  THE  1ST  ARRAY 

LAST  VALUE  OF  NIST  IN  IST-ARRAY  :  DIMENSION  * 

OF  1ST  * 

PARAMETER  TO  TELL  WHETHER  TO  STORE  FOR  TABULATION 
<IKEEP  =  1)  OR  PASS  ON  TO  NFXT  SECTION  ( IKEE  F  =  0 )  * 

************************* 
MEAN  FLOW  DEPTH  (FT)  * 

H  NONQIMENSIONALIZED  BY  H - THUS  UNITY  * 

H  "  "  W  * 

RADIUS  OF  CURVATURE  (FT)  * 

RC  IS  POSITIVE  IF  THE  ORIGIN  IS  LOCATED  ON  THF  * 
RIGHT  BANK  SIDE 

PC  IS  NEGATIVE  IF  THE  ORIGIN  IS  LOCATED  ON  THE  * 
LEFT  BANK  SIDE  * 

PREVIOUS  SECTION’S  RC  TO  CHFCK  ORIENTATION  * 

STARTING*  ENDING  CENTERLINE  COORDINATEc  OF  EACH  * 
BEND  SEGMENT  * 

CUMULATIVE  CENTERLINE  STREAMWISE  COORDINATE  * 

PREVIOUS  SECTION’S  CENTER  LINE  COORDINATE  * 

CENTER  LINE  DISTANCE  BETWEEN  SECTIONS  * 

DARCY-WE ISBACH  FRICTION  FACTOR  * 

NOTE  THAT  F  IS  USUALLY  DETERMINED  BY  « 

F  =  8.0  *  G  *  R  *  S  /  (  V**P  )  * 

THUS*  SPECIFYING  3  OF  (  F*  R*  S*  OR  V  )  * 

DETERMINES  THE  FOURTH  PARAMETER  * 

MEAN  STREAMWISF  FLOW  VELOCITY  (FT/SEC) 

REMAINS  IN  ORIGINAL  UNITS  BY  THE  NAME  VACT  * 

CENTERLINE  WATER-SURFACE  SLOPE  * 

RIVER  WIDTH  (FT)  * 


B2 


REMAIN'S  IN  THE  ORIGINAL  UNITS  BY  THE  NAME  UACT 
N  =  NUMBER  OF  STREAM  WISE  POSITIONS 

M  =  NUMBER  OF  RADIAL  POSITIONS - TAKEN  ODD  FOR 

CENTERLINE 

NS  L  G  =  NUMBER  OF  CONSECUTIVE  SEGMLNTEO  BENDS 
NSTFPS  =  NUMBER  OF  SECTIONS  THAT  INTERVAL  (  SI*  S2  )  IS 
TO  3E  DIVIDED 

««»**•*  ***********  ******  ***** 

ALPHA  =  CONSTANT  USED  IN  BED-LAYER  RELATIONSHIP 
BETA  =  CONSTANT  USED  IN  SHEAR-STRESS  RELATIONSHIP 
POR  =  PE D -LAYF R  POROSITY 
SG  =  SPECIFIC  GRAVITY  OF  SEDIMENT 
D  5  C  =  MEDIAN  BED-MATERIAL  PARTICLE  SIZE  <  MM  ) 

R  MU  =  DYNAMIC  VISCOSITY 
THETAC  =  SHIELDS  PARAMETER 

***************************** 
DIMENSTON  VI  <  1 7  )  »  VIM1I17)*  VNI(17) 

DIMENSION  UPIC17)*  UBIM1(17)*  UBN I ( 1 7) *  UI(17) 

DIMENSION  SI  ( 1  7  )  «  SIMl(17>t  HLIV<17)*  HLIM1VC17) 

DIMENSION  HL  IH ( 1 7 ) *  HLIM1H(17) 

DIMENSION  R  <  17 )  *  V  A  <  1  7  )  *  VQC17) 

DIMENSION  T  UB  < 1 7 ) «  TEMV<17>*  THLV<17>*  THLHC17) 

DIMENSION  TEMPR (17) «  TEMPS(17) 

DIMENSION  I  L  OC  <  5  )  «  I  ST  C  6  ) 

DIMENSION  V (6*17)*  UB<6*17>*  U<6*17),  HL<6*17) 

DIMENSION  UDPU (6*17)*  UTRAN<6*17)*  ANGLE<6*17) 

DIMENSION  GOSS <6*17) 

P&IMOS  I/O  COMMENT 

OPENtDfFILE-’SFGDAT*) 

0PrN(6*FILE=,0UTT«) 

READ  <5, IP)  VBAR  *H*W*SCL  *POR  *SG*RMU*NSEC 
URI TE<6*b) 

5  FOR MAT(//*20X, ‘MATHEMATICAL  MODEL  FOR  THE  PREDICTION  OF  • 
t  */ *  2  4  X  « •  TH  E  VELOCITY  FIELD  IN  RIVER  FLOW**//) 

WRITE <6* 6) 

6  FORMAT  <1X*7B(*** )) 

10  F0RMAT<6F1? ,5*E15.b* 14) 

NSFSP1  =  NS EG  ♦  1 

READ  (5*12)  <  ILOC(I  )*I=1*NSEGP1  ) 

12  F0RMAT<7I1J> 

WRITE(fe»14) 

14  FOR  «AT< //*1BX,»VALUES  FOP  VBAR*  H*  U*  SCL*  POR*  SG  *  R  ML)  * 

J  NS EG  :•) 

WRITE<6*15>  VBAP,  H*  U*  SCL*  POR*  SG*  R  MU*  NS  EG 
lr-  FOPMAT(5X,3F10.3,E13.5*2F7.i»E11.3,I4*/) 

WRI TE<6,16) 

16  FOR  MAT  <  14X*  • SECTION'  NUMBERS  WHF»E  MEW  SEGMENTS  BEGIN 
t  auE  : • ) 

WRI TE(6,12)  <  ILOC< I ) *  I = 1  * NSE GP 1  ) 


noomnnnonnonnnnnnnorjniionnnnnorinn 


WRITECf *6) 

C*«** ****«*,*+*•*******»••*•, 


DEFINITIONS  OF  VARIOUS  TERMS  USED 


RADIAL  CENTERLINE  POSITION 

DESIGNATED  RADIAL  POSITION  TO  RE  PRINTED 


MP  = 
NP  = 
I  = 
TOUT  = 

ROUT  = 
EPSV  = 
EPSU  = 
K  CUNT  = 
KM  AX  = 

KTIM  = 

K°  = 
KCPT1  - 


KCPT2  = 


K  OP  T  3  = 


NOTE:  5  POSITIONS  PRINTED  IN  ALL:  * 

1  =  NEGATIVE  RIGHT  BANK  * 

?  =  QUARTER  POINT  * 

3  =  CENTER  LINE  * 

4  =  3-QUARTER  POINT  * 

5  =  POSITIVE  LEFT  BANK  * 

INCREMENT  FOR  RADIAL  LOCATION  OUTPUT  * 

INCREMENT  FOR  DOWNSTREAM  SECTION  OUTPUT  * 

DOWNSTREAM  SECTION  INDEX  * 

D/S  PRINTING  OPTION  PARAMETER  FOR  RESULTS  * 

EVERY  IOUT-TH  SECTION  IS  PRINTED 

TEST  PARAMETER/VARIABLE  FOR  IGUT-TH  SECTION  * 

EPSILON  FOR  ERROR  IN  V  BETWEEN  ITERATIONS  * 

*  "  "  •  UBAR  BETWEEN  ITERATIONS  * 

ITERATION  COUNTER  OF  UBAR  * 

MAXIMUM  NUMBER  OF  ITERATIONS  DESIRED  FCR  UBAR  * 

CALCULATION 

FREQUENCY  OF  OUTPUT  DURING  ITERATION  PROCEDURr  FOR* 
UBAR  * 

PRINTING  PARAMETER  FOR  KTIM-TH  ITERATION  OUTPUT  * 
OPTION  PARAMETER  FOR  INITIAL  APPROXIMATION  FOR  UBAR 
1  FOR  USING  PREVIOUS  SECTION'S  VALUES  * 

0  FOR  USING  DISCRETIZED  CONTINUITY  EQUATION  * 

3  FOR  ANALYTIC  SOLUTION  OF  THE  CONTINUITY  EQUATION* 
OPTION  PARAMETER  FOR  SOLVING  QUADRATIC  FORMULA  * 
FOR  V  * 

1  FOR  ORIGINAL  DISCRETIZED  MOMENTUM  EQUATION  * 

*>  FOR  MODIFIED  DISCRETIZED  MOMENTUM  EQUATION  WHICH* 
INCLUDES  THE  CONTINUITY  EQUATION  * 

OPTION  PARAMETER  FOR  RADIAL  INTEGRATION  DIRECTION  * 
1  FOR  DIRECTION  FROM  THE  POSITIVE  LEFT  BANK  ACROSS 
TO  THE  NEGATIVE  RIGHT  BANK  * 

?  FOR  DIRECTION  FROM  THE  NEGATIVE  RIGHT  BANK  ACROSS 
TO  THE  POSITIVE  LEFT  BANK  * 


MORr  DEFINITIONS 

=  PI  !!!  BOSTON  CREME  *  APPLE*  PUMPKIN*  t T  C ♦ 
i  =  GRAVITATIONAL  CONSTANT 
t  =  COEFFICIENT  IN  SEDIMENT  RELATION 
<  =  EXPONENT  IN  SEDIMENT  RELATION 

SEDIMENT  RELATION  -  POWER  LAW  OF  THE  FORM 

QS  =  AAA  *  VBAR  **  BB'< 
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FT  MM 
KUVM 
ISTP 


CONVERSION  NUMBER  OF  MM  IN  ONE  FOOT  * 
MAX  NUMBER  OF  ITERATIONS  FOR  U-V  COMPATIBILITY  * 
SIZE  OF  IST-ARRAY - PR  T  NT  It  ISTD - LAST  VALLE* 


OS 
QSACT 
V  STAR 
FR  CUOE 
RU 
PEL* 
BL  *  PR 


=  WATER  DISCHARGE 
=  SEDIMENT  DISCHARGE 

=  SEDIMENT  DISCHARGE - ACTUAL 

=  SHEAR  VELOCITY 
=  FROUOE  NUMBER 
=  POWER  LAW  EXPONENT 
=  RADIAL  STEP  SIZE 

=  POSITIVE  LEFT  AND  NEGATIVE  RIGHT  BANKS* 
RESPECTIVELY 


**************.******«****** 


C  GP 1 tCP?  =  DIMENSIONLESS  GRAVITATIONAL  TFRMS  BY  DEPTH  H  AND 
C  WIDTH  U 

C  RESTAR  =  BOUNDARY  REYNOLDS  NUMRER 
C  VSTARC  =  CRITICAL  SHEAR  VELOCITY 
C  FRO  =  DENSIMETPIC  FROUDE  NUMBER 


V  AR I ABLr  S  TO  DEFINE  EVERY  RUN 
PI  =  3.1415B2654 
G  =  32.174 
rTMM  =  304.8 
AAA  =  0.108 
BBfi  =  a.O 
N  =  54D 
»  :  17 
MCL  =  T 
I  Ml  =0 
IM?  =  MCL 
IM?  =  13 
IM4  =  M 
IOUT  =  2 
MP  =  1 
NP  -  1 

F  PS  V  =  0.00  1 
F  PS  U  =  C.01 
KM A  X  r  20 
KTIM  =  10 
KOPT1  =  3 
KOPT2  =  2 
KOPT3  =  2 
KUVM  =  20 

I  ST  0  =  6 

1ST  Cl)  =  1 
1ST  <  2 )  =  60 
1ST  <  3)  =  2GB 
1ST (4)  =  341 
1ST ( 5  )  =  477 
1ST (ISTD)  =  N 


% 


WRITE<6,20>  N *  M,  MCL,  IM1,  IM2,  IM3,  IM4*  IOUT,  MF,  MR 

23  FOR MAT<5X, ’DOWNSTREAM  STEPS  =  »,I5,’  RAOIAL  STEPS  =  ’,15, 

J  •  CENTER  AT  M  =  • *  1 5 * / , 5X , * R AO  I AL  POSITIONS  c  T  OR  r  D  AT 

S  J  =  1  *, 415, /,5X, ’RESULTS  OUTPUT  EVERY  ‘fib* 

$  •  SECTIONS *,/,5X,D/S  AND  RADIAL  OUTPUT  FREQUENCY  IS  *» 

S  215**  STEPS*,/) 

WRITE<6,22)  EPS V t  EPSU,  KMAX,  KTIM,  KOD  T 1 »  KOPT?,  KOPT  1 
22  FORMAT (5V, ’RELATIVE  ERROR  CRITERIA  FOR  V  AND  UHAR  ARE  ’, 

*  2F10. 5, /,5X, ’MAXIMUM  ITERATIONS  »,I5,»  PRINTED  EACH  ’,15, 

*  ’ITERATIONS’,  /  ,5X , ’PR OGR AM  OPTIONS  FOR  URAR,  "OMENTUM 
t  FORM  * D IRECTION  ARE  •,  315,/) 

WRITECfe,24)  NN,  KUVM,  ISTD,  ( IS T ( K ) ,K  =  1 , I  ST D  ) 

24  F0RMATI5X,* INITIAL  NUMBER  OF  SUBINTERVALS  FOR  SIMPSON 

*  RULE  ♦  * , 1 5 ,  /,5X,«NAX  NUMBER  OF  U-V  ITERATIONS  IS  =  *, 

S  I5*/*5X*  ’DIMENSION  OR  NUMBER  OF  SECTIONS  TO  Pt  TABULATED 
t  IS  =  **I5,/*5X,  • AND  ARE  AT  SECTIONS:  • * < / *5 X « 7  1 1  G > > 
URITE<6,25)  AAA,  BP  B 

25  FOR " AT  < / ,5X, ’SEDIMENT  POWER  LAW  OF  THE  FORM  QS  =  A  * 

%  <  V  )  **  B  * ,  /,5X,*  WITH  A,  B  =  ’,2F12.4,/) 

WRI TE(6,6) 


C  COMPUTE  QUANTITIES  TO  BE  USED  IN  THE  PROGRAM 

Q  =  V3AR  *  H  *  W 
QSACT  =  AAA  *  VSAR**PBD 
VSTAR  r  SQR  T  <  G  *  H  *  SCI  ) 

FROUDE  =  VBAR  /  SQR  T (  G  *  H  ) 

F  =  fl.O  *  G  *  H  *  SCL  /  C  VBAR**2) 

RN  =  1,0  /  SORT  <  F  ) 

RN2  -  (  RM  ♦  1,0  )**2  /  <  RN  *  (  RN  ♦  2,0  )  ) 

DEL  R  =  W  /  <  M  -  1  ) 

BL  =  W  /  2.0 
3R  =  -W  /  2.0 

WRI TE(G,b5)  Q,F,RN, DELR ,EL,fiR, VS  TAR , FROUDE, RN2,GSA  CT 
65  F0RMAT<5X, ’DISCHARGE  =  ’»F13.4,»  D AR CY - WE  I S B ACH  F  =  ’, 
%  E14.6,/ *5X»  ’POWER-LAW  N  =  *»E16.8*’  RADIAL  STEP  =  *, 
J  E1?,4,/,5X,  ’LEFT  R  RIGHT  RANK  AT  R  =  • *2E 12 ,4 ,/ , BX, 

1  ’SHEAR  VELOCITY  *,E14.6,*  FROUDE  NO  =  ♦ , E 1 4 .6 , / , 5X, 
i  ’N-TERM  GIVEN  CY  • , E 1 4 . 6 , / ,5 X, » QS  =  *,F15.6,* 

I’ENG  TONS/DAY’,/) 

WRI TE<6»6) 


C  DETERMINE  THE  NOND IM EN SI  ON AL  TERMS 

VACT  =  VBAR 
WACT  =  U 

GPI  =  G  *  H  /  VBAR  *  *  2 
G P?  =  G  *  W  /  VBAR  * *2 


►y 

V 
.  « 


HH  =  H  /  w 
BL  =  BL  /  W 
SR  =  HR  /  u 
DELR  =  DELR  /  W 
VBAK  =  1.0 
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\  A*. 


l.C 
1.0 
1.0 

UR1 TE(b*70) 

70  FORMAT*//, 2 7 X , ’MONO IMENS I  ON AL  OUANTITIES  :  •) 

WRITE(6,72>  GP1,  GP2,  HH,  BL,  BR*  DELR 
72  FOR  M  AT ( 5X  **  TWO  GRAVITY  TERMS  FOR  -VERT-,  -HOR-  ARE  *, 

J  2E14.b,/,5X,  ’DEPTH  =  ’.E14.6,*  LEFT  AND  RIGHT  BANKS  AT  * 
S2E14.B,/,5X,  ’RADIAL  STEP  ’,E14.6> 

URITECf ,6) 


DEFINE  THE  INITIAL  ORIENTATION  OF  THE  RADIAL  POSITIONS 
I F (  KOPT3  .EQ.  2  >  GO  TO  7b 

RADTAL  POSITIONS  DEFINED  FROM  THE  POSITIVE  LEFT  BANK 
TO  THE  NEGATIVE  RIGHT  BANK 


**************  K0PT3  =  1  OPTION  ************************* 

DO  74  d  =  1,  M 

R  < J )  =  BL  -  <  J  -  1  >  *  DELR 

74  CONTINUE 
GO  TC 

75  CONTINUE 

RADIAL  POSITIONS  ARE  DEFINED  FROM  THE  NEGATIVE 
RIGHT  BANK  TO  THE  POSITIVE  LEFT  BANK 

**«•*•«****•**  KOPT?  =  2  OPTION  ************************* 

DO  78  J  =  1,  M 

RCJ)  :  BR  ♦  (  J  -  I  )  *  DELR 
7R  CONTINUE 
70  CONTINUE 


SEG  =  1 
=  0 

IIST  s  1 

***«•«*•*  NEW  SECTION  ********************************* 

«0  CONTINUE 
1  =  1+1 

I F  <  1ST (NIST)  .EG.  I  )  GO  TO  95 
CO  NOT  WANT  TO  STORE  THIS  SECTION’S  DEPTHS  AND  VELOCITIES 
IKFFP  =  0 
GO  TO  100 
c 5  CONTINUE 

WANT  TO  STORE  THIS  SECTION’S  DEPTHS  AND  VELOCITIES  FOR  LATER 
TABULATION 


IKrEF  =  1 
IK  =  NIST 
NIST  =  NIST  ♦  1 
100  CONTINUF 

I F (  I  .EC.  2  )  GO  TO  220 
IJK  =  I L  0  C ( I SEG)  ♦  1 
I F (  I  .CQ.  IJK  )  GO  TO  120 
I F (  I  .NE.  1  )  GO  TO  220 
120  CONTINUr 


*  *  *  *  + 


********  EITHER  HIVE  SECTION  DEEP  INTO  A  NEW  S  r  G  M  E  N  T 

OR  AT  THE  INLET  SECTION 

ISEG  =  I SEG  ♦  1 
REA  5 (5,125)  RC,  SI*  S 2*  ALPHA*  BETA*  THFTAC*  PHC,  LSTEFS 
125  F0RMAT(7FU.5,I3) 

WRI TE (b, 6) 

URITE(f>,127)  I *R  C*S 1 *S2  * ALPH A  *  BE  TA *  THE T AC  *  05 0  *  NS TE PS 
127  FOR  PAT (//»17X»  •  NEW  SEGMENT-IMPORTANT  PARAMETERS  GIVEN  AS  • 

5  */*  5X,*F0R  SECTION  I  =  *,I5,*  RC  =  • ♦ El  A .6  * / * 5X * »SE bME N T 
5  LOCATION  BETWEEN  » * 2E  1 5 . 6 */ *5X *  * AL PH A  *  BETA  =  f*2F8*4, 

$  •  THETAC*  D50  =  **  2F8 .4 */ ,5X *• NUMBER  OF  SECTIONS  BETWEEN 
SSI  R  S?  IS  * , 15 ♦ // ) 

COMPUTE  OTHER  VARIABLES 
RESTAR  =  VSTAR  *  050  /  (  RMIJ  *  F TMM  ) 

VST  ARC  =  SGRT(  (SG  -  1.0)  *  G  *  D50  *  THETAC  /  r  TM  M  ) 

RATIO  =  VSTAR  /  VSTARC 

FRO  =  VACT  /  SORT (  (SG  -  1.0)  *  G  *  050  /  F TMM  ) 

OS  =  (  S2  -  SI  )  /  NSTEPS 

COMPUTE  NONDIMENSIONAL  QUANTITIES 
RC  =  RC  /  VACT 

51  =  SI  /  VACT 

52  =  S2  /  WACT 
OS  =  OS  /  U ACT 

05C  =  050  /  (  F  T MM  *  H  ) 

CALL  PG(F,VBAR,050»THETAC»POR*ALPHA*BETA,SG*CP1 *G1  ,G2*C3) 
NOTE:  G1  =  Gl<  N *  F.  BETA  > 

G2  =  32 (  N  ) 

G 3  =  G 3 (  BETA,  ALPHA,  POR,  F,  THETAC,  VBAR, 

G,  05C,  DRHOS  ) 

URI TEC6, 130) 

130  FORMATt// »5X»*C0MPUTED  VARIABLES  FOR  THE  NEW  SEGMENT 
S  GIVEN  BY*,/) 

URITE(6»132)  RESTAR,  VSTARC,  RATIO,  FRO,  Gl,  G2  ,  G3 
132  FORMAT  <5X,*REST  AR,  VSTARC*  RATIO  =  •  ,  3E  1 5 .  G  » /  , 5 X  ,  *  PF.  NS  T  - 
SMETIC  FROUOE  =  **E14.6,/,5X»*G1,  G? ,  G.S  =  »,3E15.7,/) 

URI TE(6,6) 

WRTTE(f , 135) 

135  F0RMAT(5X,*N0N0IMENSI0NALIZED  QUANTITIES  GIVEN  PY  •*//) 

URI TE (6*137)  RC,  SI,  SP,  D50*  OS 

137  F0RMAT(5X**RC,  SI,  S?  =  *,3E14.6,*  D50  *  * E 1 4 . 6 , / , 5 X, 
t  ’INTERVAL  BETWEEN  SECTIONS  OR  =  *,E15.7,/) 

WRI TC(6,6> 

TEST  FOR  INLET 

TEST  FOR  REQUIRED  NEW  RADIAL  ORIENTATION 
I F (  I  .EG*  1  )  GO  TO  145 
OR  I :  N  =  RC  /  R J 
R  0  “  PC 

I F  f  ORIEN  .GE.  0.0  )  On  TO  220 

138  CONTINUE 

I F (  KO? T 3  *EQ.  1  )  GO  TO  1 3 3 


GO  TO  140 
1  39  CONTINUE 
KOp  T3  =  2 
140  CONTINUE 

00  14?  J  =  It  M 
TEV  °R ( J  )  =  R  <  J ) 

TFMPSCJ)  =  SIM1 ( J) 

TUOCJ)  =  UR  I  Ml (J) 

TEVV(J)  =  VIMl(J) 

THLV(J)  =  HLIMIV(J) 

THLH(J)  =  HLIMIHIJ) 

142  CONTINUE 

DO  144  J  =  It  M 
JJ  =  <  M  ♦  1  >  -  J 
R ( J )  =  TEMPR(JJ) 

SIMI(J)  =  TFMPS(JJ) 

VIMl(J)  =  TEMVfJJ) 

U3IMKJ)  =  TUB  <  J  J ) 

HLIMIV(J)  =  THL V ( JJ ) 

HLIMIH(J)  =  THLH(dJ) 

144  CONTINUE 
GO  TO  220 

145  CONTINUE 

C  ************************  INLET  SECTION  ******************* 

C  SET  ARBITRARY  CONDITIONS - CAN  IMPOSE  ANYTHING 

R  0  -  R  C 
SC  =  0.0 
UO  =  0,0 
160  CONTINUE 

C  ESTIMATE  INITIAL  DOWNSTREAM  V-VELOCITIES  AT  INLET 

C  SECTION  USING  DARCY-WE ISBACH  APPR OX  I M AT  I  ON - NO TE  THAT 

c  VALUES  at  BANKS  ARE  NONZERO 

DO  170  J  =  1*  M 
HLIV(J)  =  HV 
HLIHCJ)  =  HH 
SICJ)  =  0.0 
UI( 0)  =  0.0 

VNI <  J  >  =  SQRT<8.0*GP1*HLIV(J>*SCL*RC/<F*(RC  ♦  R<J>>>> 

170  CONTINUE 

WRITrift* 174) 

174  FOR  MAT ( /  *  5X  t  ’INLET  SECTION  V  =  •*/) 

WRITF<6»176>  tVNKJ)*J  =  l*M*MP) 

176  F0RMAT(/*(3X«5E15.7>) 

«RI TEC6* 177) 

177  FORMAT<//) 

r  CALCULATE  THE  FLOW  DISCHARGE  AND  THEN  ADJUST  THF  V- 

C  VELOCITIES  TO  SATISFY  THE  MEAN  FLOW  CONTINUITY 

CALL  PON(MfR,VNI*HLIV*VA,VO*AT,OT) 

F  TA  =  1.0 
ETA  =  (3  /  QT 
WRTTE<6*178)  Q«  OT*  ETA 


178  FORMAT  ( //,5X,*Q,  OT,  ETA  =  * ♦ 3E 1 6 . 8  * // ) 

DO  180  J  =  1  ♦  M 
VNI(J)  =  VNHJ)  *  ETA 
180  CONTINUE 

CALL  PQNCM,R,VNI*HLIV«VA,V(UAT,QT> 

WRITE<6, 182) 

182  FORMAT  <5X** ETA-MODIFIED  V-VELOCITIES  WITH  UMIT 
J  01 SCHARGES • ) 

DO  ISO  J  =  1,  M,  MP 

WRITE (6*185)  J  »  VN I (J)«VA(J)«VG(J) 

185  FOR M  AT  (5X  *  *  J  =*,I3t*  V  =  **C16.8**  VA  =*,E16.8,*  V  (1  =  ♦, 

T  E16.8) 

ISO  CONTINUE 

WRITE(Ff 1^1)  QT 

191  FORMAT (5X  ,  *NEW  VALUE  OF  QT  WITH  MODIFIED  V  -  **E16.8> 

C  COMPUTE  SEDIMENT  DISCHARGE 

QS  -  0.0 
DO  193  J  =  If  M 
VVV  =  VNI(J)  *  VACT 
RATA  r  V  A ( J )  /  VA(MCL) 

QS  =  QS  ♦  RATA  *  VVV**BBB 

193  CONTINUE 

QS  =  QS  *  AAA  /  <  M  -  1  ) 

QSNO  =  QS  /  QSACT 
WRITE (6,194)  QSt  QS ND 

194  F0RVAT(/»5X, *QS  s  *,E14.6f»  ENG  TONS/DAY  QS/QSACT  =  », 

J  E15.6,/) 

C******************************* 

C  FIND  INITIAL  VALUES  FOR  UBAR  AT  INLET  SECTION 

C  FROM  THE  ANALYTICAL  SOLUTION  OF  THE  SIMPLIFIED 

C  CONTINUITY  EQUATION  USING  A  ZERO-VALUE  BOUNDARY 

C  CONTINUITY  AT  THE  OUTSIDE  BANK 

C  USE  CONTINUITY  EQUATION  WITH  D AR C Y-W E I  SB ACH  FOR 

C  0  <  V  *HL ) / DS  AT  INLET  TO  APPROXIMATE  UBAR 

TEMP  =  O.C 
UONIQ)  =  0.0 

IF ( KCPT  3  .EQ.  1)  TJM1  =  (-2.0+BL/RC)  *  SORT ( 1 . 0 ♦PL /RC > 

IF(K0PT3  .EQ.  2)  TJM1  r  (-2.0  +  BR/RC)  *  SQR T ( 1 . 0 ♦ BR /RC ) 

Cl  =  8.0  *  GP2  *  SCL  /  (  P  *  MH  ) 

Cl  =  -1.0  ♦  ETA  *  G2  *  G3  *  SORT  (  Cl  )  *  Rf**2 
DO  200  J  =  2,  M 
RADR  =  RC  ♦  R ( J ) 

UAL  =  <  -2.0  ♦  R ( J )  /  RC  )  *  SORT (  1.0  ♦  R(J)  /  RC  ) 

UBNI(J)  =  (  TEMP  ♦  Cl  *  (  UAL  -  TUMI  )  )  /  RADR 

TEMP  =  U  BN  I  ( J)  *  RADR 
TJMi  =  UAL 
200  CONTINUE 

URITE(6f210)  I,  (U'JNI  <U)  ,J  =  1  ,m,MP) 

210  FORwAT(/f»  I  =  *  ,  1 3  ♦  *  VALUES  FOR  UBAR  =  • * / , ( 4 X . 5E 1 c . 7 ) ) 
WRI TE<6»  212 ) 

2 12  FOR MAT<//»5Xf *ENO  OF  INLrT  SECTION**//) 


BIO 


WRI  TE<6,6) 
GO  TO  763 
22U  CONTINUF 


C  *******************  NEW  SECTION  ******************* 

C  SAME  O'?  NEW  SEGMENT 

C  FOR  THE  NEW  SECTION*  DEFINE  THE  ANGULAR  COORDINATE 

C  AND  THE  STREAMWISE  POSITIONS  ACROSS  THE  TRANSVERSE* 

C  SI.  THEN*  DETERMINE  THE  CENTERLINE  SECONDARY  FLOW 

C  VELOCITY  UCLI*  AND  THE  TRANSVERSE  BED  SLOPE*  ST1* 

USING  THE  CENTERLINE  STREAMWISE  DISTANCE. 

FINALLY*  CALCULATE  THE  LOCAL  FLOW  DEPTH  AT  EACH  OF 
THE  RADIAL  POSITIONS. 


2  21  CONTINUE 
sen  =  SC 
SC  =  SC  ♦  OS 

UU1  =  UO  /  E  XP  <  G1  *  C  SC  -  sen  )  /  HH  ) 

UU2  =  G2  *  VLAR  /  EXPC  G1  *  SC  /  HH  ) 

UU3  =  HH  *  (  EXP ( G 1 *SC/HH>  -  E X P < G 1 *SC 0 /HH )  >/C  G1  *  RC  ) 
UC  =  UU1  ♦  UUP  *  UU 3 
STI  =  UC  *  G3  /  VBAR 
UO  =  UC 

DO  225  J  =  1*  M 

HLI H(J)  =  HH  ♦  R (J)  *  STI 

HLIV(J)  =  HLIH(J)  /  HH 

TEST  FOR  NON-POSITIVE  DEPTH.  IF  SO,  MUST  MODIFY 
PROGRAM. 

FXIT  IF  ENCOUNTERED. 

IF(HLIVCJ)  .LE.  0.0)  WSITEC6,22R>  I,  J,  HLIVCJ) 

22A  F0RMAT(//*5X*’ NON-POSITIVE  FLOW  DEPTH  OCCURS  AT  SECTIOA 
t  I  =  * , I  5  *  / *5X , ’RADIAL  POSITION  J  =  * « I  5  * 5X * • F L OW 
*  DEPTH=  *,E2.4*/*5X*  ’EXIT  FROM  PROGRAM  TO  MODIFY*) 
IF(HLIVCJ)  .LE.  0.0)  GO  TO  950 

225  CONTINUE 

DETERMINE  INITIAL  VALUES  FOR  SECONDARY  FLOW  USING 
FALCON’S  RELATION  AMD  THE  PREVIOUS  SECTION’S  V-VALUES 
DO  226  J  =  1,  M 

FAC  =  HLIV< J)*VIM1C J)*RC  /  <  HH * V I M 1 ( MC L ) * < R C  ♦  R(J>>  ) 
UI(J)  =  UC  *  FAC 

226  CONTINUE 

DANCLE  =  <  SC  -  SCO  )  /  RC 
DO  230  J  =  1*  M 

SI(J)  =  <  RC  ♦  R ( J )  )  *  DANGLE  ♦  SIMl(J) 

231  CONTINUE 

f.CTF :  FOR  RC  <  0,  THE  TWO  NEGATIVE  OUANTITITS  WILL 

STILL  YIELD  A  POSITIVE  STREAMWISE  COORDINATE*  S. 
I F ( K  OPT  1  .NT.  l)  GO  TO  235 


*  *********KnPTl=lCPTION«  *******  **«* 

FOR  A  ROUGH  APPROXIMATION,  ASSUME  THAT  THE  CURRENT 
VALUES  OF  UPAR  ARF  THE  S  AM  F  AS  THE  VALUES  CALCULATED 
4T  THE  P  °  LV I  0  US  SECTION. 


DO  232  J  =  2  *  M 
UBI  <  J)  =  UBIMKJ) 

232  CONTINUE 
GO  TO  250 
235  CONTINUE 

IF<  KGPT1  .NE.  2)  GO  TO  241 

************  KOPT1  =  2  OPTION  *********** 

AS  A  ANOTHER  ROUGH  APPROXIMATION*  DISCRETIZE  THE 
CONTINUITY  EQUATION  TG  APPROXIMATE  THE  VALUES  F  CD  U;'AR 
IN  THE  FIRST  ITERATION. 

TEMP  =  0.0 

B 1  -  2.0  *  GP2  *  SCL  /  F 

B2  =  -1.0  *  G1  *  SI(MCL)  /  HH 

B 1  -  -3.0  *  G2  *  G3  *  SDR  T  (  B1  )  *  E XP (  B2  > 

DO  240  J  =  2*  M 
DELS  =  R (d)  -  R ( J-l  ) 

TJ  =  HL I H ( J )  *  (  RC  ♦  R  < J)  > 

F  AC  J  =  R  ( J )  /  (  (  RC  ♦  RtJ)  )  *  SORT  ( T  J  /  RC)  ) 

UAL  =  R1  *  DELR  *  FACJ 
UBI ( J>  :  UAL  ♦  TEMP  /  TJ 
TEMP  =  UBI ( J  >  *  TJ 

240  CONTINUE 
GO  TO  250 

241  CONTINUE 

+  ********ft*ft  +  ft**  +  ftftftftftftftftft  *  ft  *  ft  ft 

*  *  *  *  *  *  *  »  *  »  *  *  KOPT1  =  3  OPTION  *  *  *  *  *  *  ***  «* 

ANALYTICAL  SOLUTION  OF  THE  CONTINUITY  EQUATION  TP  pr 
USED  TO  SOLVE  FOR  UBAR  AS  AN  INITIAL  APPROXIMATION  '.HEN 
2TAPTING  A  NEW  SECTION. 

Z1  =  -3.0  *  G2  *  G3  *  SOR T  <  2.0  *  GP?  *  SCL  /  E  ) 

Z1  =  Z1  *  E  X  P (  - G1  *  SI(MCL)  /  HH  ) 

TENT  r  0.0 
DO  ?4F,  J  =  2*  M 

CALL  EVAL(STI*HH*RC,RCJ>tR<J-l)»PIA) 

90TJ  =  HLIH(J)  *  (  RC  ♦  R(J>  ) 

UBI <  J )  =  <  Z1  *  RIA  ♦  TENT)  /  BOTJ 
TENT  =  U?I(J)  *  POTJ 
246  CONTINUE 
250  CONTINUF 


255  CONTINUE 

FIND  STREAMUISE  VELOCITY  V  BY  SOLVING  MOMENTUM  EQUATION 
THAT  HAS  BEEN  DISCRETIZED  INTO  A  QUADRATIC  FORM  USING  A 
SIMPLE  BACKWARD  DIFFERENCING 

READY  TO  BEGIN  SOLUTION  OF  DISCRPTI/fP  MOMENTUM  r JU A T 1 0* 
R  OF  V*  EiUT  FIRST,  NEED  A  BOUNDARY  CONDITION  FOR  V.  USE 
DARCY-WFIfBACH  EQUATION  AT  BANK. 

I F ( K  CP  T  3  .EC.  1)  RJ  =  BL 
I  F  (  K  C  P  T  3  .EG.  ?)  RJ  -  i-,~ 

VI(I)  =  SQR T (8. U*GP 1 *HL I  V (  1  > *SCL  ♦RC  /  (F*(RC  ♦  RJ))) 

VNI  <  1 )  =  VI  <  1  ) 


T  F.  W  <  1  )  =  V  1  (  1  ) 

K  0 1 J  U  T  =  n 
ETA  =  1.0 
K'JVT  =  0 
KGT  =  o 
KTT  =  0 
21?  CONTINUE 
KO  =  0 

3  If,  CONTINUE 
KUV  =  C 
220  CONTINUE 

•***«**«•*•***  DISCRETIZED  MOMENTUM  EQUATION  FOR  V  ***** 

- CENTER  HERE  FOR  NEW  ITERATION  AT  SAME  SECTION  - 

ERR  V  r  0 . 0 
DO  EUO  J  =  ?♦  M 

DETERMINE  F-FUNCTIONS  AT  SPECIFIC  POSITION  ( S  *  R  ) 

CALL  PF(F*HKfRC*VBAP*Gl*G2*G3fR(J)tSTI*UICJ)*Fl*F2*F3*F-») 

DEL S  =  S I <  J )  -  S  I M 1 (  J) 

STEPR  =  R  <  J )  -  R(d-l) 

CALCULATE  COEFFICIENTS  FDR  QUADRATIC  FUNCTION  IN  U: 

AA  *  U  *  *  2  ♦  ri  B  *  V  ♦  CC  =  0.0 
T R  s  HH  *  F  A  /  STEPP 

TS  =  HH  *  FA  /  DELS 

******************************* 

I F ( K  OPT 2  .NE .  1 )  GO  TO  340 

*  ***  *  *  *  *  *****  *K0PT2=1  *******  *  *  *  ** 

COEFFICIENTS  FOR  REGULAR  DISCRETIZED  STREAMWISE  MOMENTUM  EQ. 
DEPEND  ONI  V<I-1.J>*VNEWCI*J-1>.U(I.J>*UBAR(I,J),UBAR<I,J-1> 
AA  =  RN?  *  (  F?  ♦  TS  >  ♦  F  /  8.3 

BB  =  (FI  ♦  TR)  *  (UF'ICJ)  ♦  UI(J)  /  (2.0  *  P  N  ♦  1)  ) 

CC  =  GP1  *  HV  *  SCL  *  FA  *  RC  /  (RC  ♦  R(J)) 

CC  =  CC  ♦  K N 2  *  TS  *  VIMl(J)**? 

Cv  r  TR  *  (  UBHJ-l)  ♦  UI(J)  /  (  2.3  *  RN  +  1.0  )  ) 

CC  =  -1.0  *  (CC  ♦  CA  *  VNKJ-l)) 

GO  to  360 

******************************** 

3  AO  CONTINUE 

***************  KOPT2  =2*********** 

COEFFICIENTS  FOR  MODIFIED  DISCRETIZED  STREAMWISE  MOMENTUM 
FQUATION  WHICH  INCLUDES  THF  CONTINUITY  EQUATION 


NOTE  : 

NO  RADIAL  DIFFERENTIATION 

IN 

UBAR 

» 

1 

1 

E  p  E  N  D 

on:  V  (  I 

VNF  U  (  I 

»  j-i  > 

» 

U( I,J) 

* 

UDAR ( I, J) 

A  A  *  F  2  * 

(  R  N  2  - 

1.0  )  ♦  TS  * 

(  RN2 

- 

0.5  ) 

♦ 

F  /8 .0 

PR  =  TR  * 

(  JBI(J) 

*  UI  ( J )  /  ( 

2.0  * 

RN 

♦  1. 

0 

)  ) 

R  B  -  p  P  ♦ 

U  I  (  J )  * 

El  /  (  2.0  * 

R  N  ♦ 

1  . 

0 

> 

CC1  =  DPI 

*  HV  ♦  SCL  *  FA  *  RC 

/  (  RC 

♦ 

R  (  J) 

) 

CCD  =  T  R  * ( 

UBI ( J  > 

♦  U I  ( J  )  /  (  2 

•  3  *  R  N 

♦ 

1.0  ) 

) 

*  V  N I  (J-l  ) 

CCA  =  T  S  * 

(  R  N  2  - 

0.5  )  *  V I m 1 ( J ) *  * 

? 

CC  =  -1.0 

*  (  CC1 

♦  CC2  ♦  CCi  ) 

***** 

*  *  *  * 

******* 

*  ♦ 

* 

* 

*  * 

* 

*  *  *  *  * 

260  CONTINUE 


CHECK  =  B8*»2  -  4.0  *  A  A  *  CC 
IF (CHECK  .LT.  G .  '> )  GO  TQ  400 

VNI(J)  =  (  -:3P  *  SORT(CHECK)  )  /  (  2.0  *  AA  ) 

C  SUM  THE  DIFFERENCES  BETWEEN  THE  OLD  AND  NEW  VALUES--- 

C  NOTE  THAT  THE  BOUNDARY  VALUE  IS  NOT  INCLUDED  IN 

C  THE  ERROR  SINCE  IT  IS  F I  YE  D 

FRRV  =  ERR V  ♦  ADS (  VNI(J>  -  VI(J)  ) 

I F  (  VNICJ)  .GE.  0.0  )  GC  TC  5  00 

c  ************  NEGATIVE  streamwise  velocity  ************** 
URITE<6.393)  It  J.  VNI(J) 

2  9  0  P0RFAT</*5X»  ’NEGATIVE  STREAMWISE  VELOCITY*./.?*. 

4  ‘SECTION  I  =•♦  15.*  RADIAL  STEP  J  =  *»T5»‘  V  =  •♦"15.7, /> 
GO  TO  415 
400  CONTINUE 

URITE(6,410)  CHr  CK 

410  FORMAT (// .5X. ‘NEGATIVE  RADICAL  IN  U-OUADRATIC - CHECK  r  *. 

4  El  0  •?./  ) 

415  CONTINUE 

C  EXIT  FROM  PROGRAM  FOR  A  NEGATIVE  RADICAL 

GO  TO  950 
500  CONTINUE 

C  ********  find  NEW  U-VALUES  BASED  ON  V  JUST  CALCULATED  ***** 

DO  505  J  =  1,  H 

FAC  =  HLIV( J)*VNI(J)*RC  /  ( H V ♦ VN I < N C L  )  *  (  RC  ♦  R(J>>  ) 

U I  (  J  )  =■  UC  ♦  FAC 
505  CONTINUE 

ERR  V 2  =  (  FRRV  /  <M  -  1)  )  /  VD  A  R 
KUVT  =  KUVT  ♦  1 

I F  €  ERR  V 2  .LE.  EPSV  )  GO  TO  510 

C  ****»•«**************•****«*••**•*****»*******»•»•.»•»..**.». 

C  RESULT  OF  V  CHANGING  TOO  MUCH  -  TRY  NEW  U-VALUES  IN 

C  CALCULATION  OF  N F VI  V-VALUES  UNTIL  CHANGES  ARE  MINOR 

K(JV  =  KUV  ♦  1 

C  REPLACE  OLD  VALUES  WITH  NEW  VALUES  AND  SOLVr  FOR  V 

C  AGAIN  BUT  NOW  WITH  THE  NEWLY  UPDATED  V-VALU^S.  NOTE 

C  THAT  U-80UNDAR Y  REMAINS  UNCHANGED. 

DO  508  J  =  2.  M 
VI ( J)  =  VN I ( J) 

V  N I ( J>  =  0.0 
5  CP  CONTINUE 

I F  f  K UV  .LT.  KUVM  )  GO  TO  32n 
WR1 TE(6*509)  I,  KUV.  EE  8  V 

509  F0RMAT(/*5X, ‘MAXIMUM  NUMBER  OF  ITERATIONS  cOR  L-V  EXCEEDED* 
I  *  5  X ,  /,5X. ‘SECTION  I  =  *.15.*  KUV  =  ‘,I5.‘  FRRV  =  * 
<E15.6./) 

GO  TO  950 
511  CONTINUE 

C  *******************  v  -  U  COMPATIBILITY  ***************** 

CALL  PON(Mf  R.VNI  .HLIV.VA.Vi.AT.i'T) 

PATO  =  A  B  S  <  0  -  OT  )  /  0 
KOT  =  KQT  ♦  1 
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ETA"  =  lTA 
ETA  =  0  /  AT 
DO  5  35  J  =  2* 

TEMV(J)  =  V\'1<J>  *  ETA 

5  35  CONTINUE 

CALL  PON(M,R*TEMV,HLIV*VA*VP*AT,GT) 

C  *•«»**..*•  FIND  NEW  UDAR  FROM  ACTUAL  IV-H)  VALUES  ****** 

ERR  u  =  0.2 
T  f  ‘1 1"  =  0.0 
U  F  N  I  < 1 )  =  0.0 
KOUMT  =  KOUNT  ♦  1 

I F ( KOPT  3  .EQ.  1)  T1  =  (  RC  ♦  PL  >**2 

IF(K0PT3  .EQ.  ?)  T 1  =  (  RC  ♦  PR  )**2 

C  USE  CONTINUITY  EQUATION  WITH  ACTUALLY  CALCULATEC 

C  VALUES  FOR  0<V*HL)/DS  AT  THIS  SAME  SECTION 

C  TO  GET  THE  SECOND  APPROXIMATION  OF  U6AR* 

C  AGAIN  WITH  UB  AR  EQUALS  ZERO  AT  THE  BANK. 

0  0  F  0  0  J  =  2  «  « 

CAPX  =  -0.5  *  <  TEMVCJ)  *  HLIH(J)  -  VIMKJ)  *  HLIMIH(J)  ) 

CAP*  =  CAPX  /  (  S I ( J )  -  SIMl(J)  ) 

T2  =  (  RC  ♦  R(.J)  )  *  *  2 

DOT  U  =  HLIH(J)  *  (  RC  ♦  R(J>  ) 

URMKJ)  =  (  CAPX  *  (  T2  -  T1  )  ♦  TFMP  >  /  POTJ 
T 1  =  T2 

TEVP  r  UF.NKJ)  *  POTJ 

ERR  U  =  ERRIJ  ♦  A  P  S  (  UBNI(J)  ■  UHIJI  ) 

600  CONTINUE 

S***.  *************.  ••***•****•»* 

C  CONVERGENCE  CRITEPIA  FOR  ERRORS:  ITERATE  UNTIL  THE  AVERAGE 

C  CIFFEREMCE  DETUEEN  ?  CONSECUTIVE  ITERATIONS  FOR  ALL  THE 

C  RADIAL  STEPS  ACROSS  THE  TRANSVERSE*  DIVIDED  BY  THE  "LAN 

C  DOWNSTREAM  FLOW  VELOCITY*  IS  SMALLER  THAN  SOME  PRESCRIBED 

C  SMALL  VALUE. 

1“************  ******  ***+».*.***•* 

ERRU 2  =  E  R  k  U  /  <  (  M  -  l  )  *  ABSC  UDNIIMCL)  )  ) 

KTT  =  KTT  ♦  l 

IF  <  £ R R U 2  .LE.  EPSU)  GO  TO  750 
IF ( K CUNT  .GC.  K M A X )  GO  TC  770 

C  *  *  *  *  *»*****»***************#****»***»***■*»*»***■»***»*******'****** 

C  REACHED  AS  A  RESULT  OF  LARGE  ERROR  WITHOUT  EXCEEDINC 

r  THE  MAXIMUN  ALLOWABLE  NUMBER  OF  ITERATIONS. 

C  NO  CONVERGENCE  YET*  SC 

C  ITERATE  THE  WHOLF  PROCESS  AGAIN  FOR  A  NEW  V  AND  UB  A  R • 

KP  =  KOUNT  /  K  T I M 
KP  =  KP  *  KTIM 
I  F ( K D  .NE.  KOUNT)  GO  TO  7  A  R 

WRI TO Cfc* 700 )  I,  KOUNT*  ETA*  ETAj*  FRRV2*  ERRL2 
7  f  i)  FOR  ”  AT  (  /  ,5X  *  *FOR  SECTION  I  =  **I5**  KOUNT  =  **I5*/,5X* 

1  •  CURRENT  AND  PREVIOUS  ETA  ARE  ♦ *2F 15 . 6  * / » 5 X* 

♦  fr?rLATIVE  ERRORS  IN  V  AND  UB  A  R  A  R  F  **2E15.6*/> 


non 


‘j 


C  WRITE  <6*  705) 

705  FORMAT(5X**PREVIOUS  VALUES  FOD  V  GIVEN  BY  •*/> 

C  WRITE! 6 *710)  !VI!J>*J=1*M> 

7 1 C  FOR  M  AT ! 5  X  *6  E 12.4 ) 

WRI TE!6*712> 

712  FOR  MAT !5X**UN-ET A -MODIFIED  MOMENTUM  VALUES  FOR 
f  V  GIVEN  BY  •*/) 

WRITECf*710)  1 VNI! J)«J=1,W> 

WRITE!6*714> 

7 1 A  FOR  ^  AT  ! 5X  * ♦ 0 LD  VALUES  FOR  UPAR  GIVEN  BY  ♦»/> 
WRITE(6»7in)  !UBI ! J) • J=1 *M> 

WRITE!6*716> 

716  F  OR  M  AT !  5X  *  *  NEW  VALUES  FOR  UBAR  GIVEN  BY  •*/> 

WRITE!  6 *710)  < UPN I ! J ) « J  =  1 1 M ) 

WRTTE!6*71R> 

718  F0RMAT(5X,*LATEST  VALUES  FOR  U  ARE  •*/> 

WRI TE!6*710)  !UI ! J) *J-1 *M) 

748  CONTINUE 

RESET  NEW  VALUES  BEFORE  REPEATING  THE  ITERATIO 
DO  740  J  =  2*  M 
VI(J)  =  V  N I ( J ) 

UBI (J)  =  UON I  ( J ) 

V N I <  J )  =  0.0 

UBNI(J)  =  0.0 

749  CONTINUE 
GO  TO  312 

750  CONTINUE 


REACHED  AS  A  RESULT  OF  A  SATISFACTORY  SMALL  ERR ' - R 
IN  UBAR  OR  V  THIS  IS  WANT  UE  WANT. 

DETERMINE  FINAL  ETA-MODIFIED  V-VALUE 
DO  752  J  =  2*  M 
VNIIJ)  =  TEMVCJ) 

752  CONTINUE 

«*<,*«•****«***•**  DETERMINE  SEDIMENT  DISCHARGE  ******* 
QS  =  0.0 
DO  753  J  =  1*  M 
VVV  =  VNI(J)  *  VACT 
RAT  A  =  V  A ( J  )  /  VACMCL) 

QS  =  OS  ♦  RATA  *  VVV** EBP 

753  CONTINUE 

QS  r  AAA  *  OS  /  <  M  -  1  ) 


NOUT  -  III  /  I  OUT 

NOUT  =  NOUT  *  I  OUT 

IFCNCIJT  »Nr  •  III)  GO  TO  763 

WRITE  (6*  755  )  I*  STI,  UC*  ERRV2*  ERRU2*  MOUNT «  KATO 
FOR  MAT (/*5X « • F  O  R  SECTION  I  =  **I4*»  ST  =  **F14.6. 

I  *U  C  L  =  **E 1 4.6  *  /  *  5  X  *  *  S  A  T  I  SF  A C  T  0 R  Y  ITERATION** 


i  GX,*ERRV2  =  *,E14.6,*  FRRU2  =  •,E14.6,/,C>X, 

*  ‘NUMBER  OF  ITfRATIONS  K  PUNT  =  *,I4,/*6X, 

i  ‘relative  difference  of  qt  to  o  ic  *,ri5.7*/> 

URI  TF  (6,757)  SC,  RC,  QT,  F  T  A »  QS«  uiSND 
767  F0KMAT(5X, ’CENTERLINE  POSITION  S  =  *,£14.6«*  WITH  R C  =  *, 
14.6,  /,5X, ’NEWEST  Q  =  *,E16.R,*  WITH  C  T  A  =  •  ,  F  1  5.  7,/ ,  CX, 
t  •  S l.  niHENT  DISCHARGE  =  *,E14.5»*  QS/QSACT  =  *,E14.5,/) 
WRI  Tf  <  *  7  f»  0  >  KTT,  KQT,  KUVT 

761  FOR  MAT ( 5X  *  *  NUMBER  OF  ITERATIONS  FOR  UB,  Q,  UV  ARE  *,316) 
WRI TE(6,71b> 

WRI  TE  (6,71 0  >  (Ur>NI(d),J=l,M,MP) 

WRI TF (6,761  ) 

761  FORwAT lbX,,ETA-MOOTFIED  MOMENTUM  VALUES  FOR  V  FIVE*  BY*) 
WRI TE (6 , 71 0 )  (VNI (J),J=1,M,MP) 

WRT  TLCb,71R) 

WPITF(6,710)  (VI ( J) , J=1 ,M,MP) 

WRI TF (6, 762  ) 

76?  F07MAT(/> 

763  CONTINUE 

c  TOR  E  VALUES  FOR  V,  UBAR,  A  N  H  U  FOR  LATER  TABULATION 
I F  <  IKFEP  .NE.  1  )  60  TO  7b 7 
I F (  K0PT3  .EQ.  1)  GO  TO  765 
HO  764  J  :  1,  V 
V  (  1  K  ,  J  >  =  V  N  I  (  J  ) 

Uli (  TK,J)  =  HUNT  (  J ) 

U ( I K  , J )  r  UI  (J) 

HL(  TK,J)  =  H  L I  V ( J ) 

OORc(lK,J)=AAA*(VNI(J)*VACT)**BHB 
UBPUriKf J)SUI( J)*UPNI< J) 

UST!RT=UBPU(IK,J)**£*VNI(J)**2 
UTR  AN( I K*J) =SQRT (USGRT) 

UC= I  I=VN I ( J ) 

IF(UCRII.ro. C#>  GO  TO  1 
ANGUU=UBPU(  !K,J)/VM  (J) 

ANGLE (IK ,J) = j7.?R578*ATAN(ANGUU) 

GO  TO  764 

1  ANGi E ( IK ,J> =90. 0 

764  CONTINUE 
GO  T0  767 

765  CONTINUE 

no  766  J  =  1,  V 

JJ  =  (  M  ♦  1  )  -  J 

V( I K,JJ)  =  VNI ( J ) 

U°(  I  K , J J )  =  UBNI (J) 

U(TK,JJ)  =  UI(J) 

ML  (  I K , J J  )  =  ML  I  V ( J ) 

G OS  S ( IK , J J ) = AA  A  * ( VN I ( J ) *  V AC  T ) **B6P 
UPPIM  IK, Jd)=UI ( J) -MIEN  I ( J) 

USO^T=UBPU( IK, JJ  )  **?*VNI (J)**? 

UTS  AN(IK,JJ)=SORT (USGRT) 

UC3  U=VNI(J) 


O  r>  ri  o  o  o 


766 

767 


768 


760 


IFCUCRII.EO.O.O)  GO  TO  3 
ANGUU=UBPU< IK*JJ)/VNI« J) 

ANGLE < IK  *JJ ) =57. 2957 8* AT AN < ANGUU) 

CO  TO  76b 

ANGLE ( IK, JJ>=90.0 

CONTINUE 

CONTINUE 

DO  768  J  =  1*  m 

HLIM1VCJ)  =  HLIVCJ) 

HLIMIH(J)  =  HLIH(J) 

SIM  1  (J)  =  SI (J) 

VIM1 ( J)  =  V  N I ( J ) 

UBIMKJ)  =  U6NI  <  J  ) 

CONTINUE 
DO  769  J  =  1,  m 
VI(J)  =  0.0 
VNI ( J)  =  0.0 
TEM  V  < J)  =  0.0 
UBN  I  <  J )  =  0.0 
UI(J)  =  0.0 
U9I ( J)  =  0.0 
HLIVCJ)  =  0.0 
HLIHCJ)  =  0.0 
SICJ)  -  0.0 
CONTINUE 


ETA 

I  F  ( 

I 

1.0 

.GE. 

N  )  GO  TO 

80  0 

1*  •# 

TRIAL 

STOPS  ********•********•**.***.** 

IFC 

I 

.GE. 

3  )  GO  TO 

95  0 

GO  TO  80 
770  CONTINUE 


REACHED  AS  A  RESULT  OF  EXCEEDING  THE  ALLOUALLF  NUMB f  R 
OF  ITERATIONS 

URITE(6*78G)  1,  KOUNT*  STI*  SC,  ERKV2*  ERRU2 
780  FOR maT(/,5X, »FOR  SECTION  I  =  **I5,*  KOUNT  =  *,I5,*  STI  =  • 
S*E1 5.6,  /,5X»* CENTERLINE  POSITION  =  *, 

t  El  4.6,/, 5X, ‘RELATIVE  ERRORS  IN  V  AND  UBAR  ARE  •,2''15.f> 
URI TE(6*716> 

WRITE(6*710)  <UBNI(J),Jrl*M> 

URI TEC6* 712 ) 

URITE(6*710)  < VNIC J), J=1 ,M> 

URI TE<6,71P> 

URITE<6*710>  CUI CJ> ,J=1 ,M> 

GO  TO  950 
POO  CONTINUE 


C  -  PRINT  FINAL  RESULTS  OBTAINED  - 

C  RESULTS  ARE  TABULATED  FROM  THE  NEGATIVE  RIGHT  BANK  TC 

C  THE  POSITIVE  RIGHT  UANK---K0PT3  =  0 
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I F (  KOPT  3  ■  E  Q  •  2  )  GC  TO  80ft 

no  ^02  j  =  l*  m 

T  EM  PR  <  J  >  =  R(J) 

P  02  CONTINUE 

DO  804  J  =  1*  M 
JJ  =  <  M  ♦  1  )  -  J 

9 < J  )  =  TEMPR  (JJ> 
e04  CONTINUE 
808  CONTINUE 


c  ESULTS  FOR  TRANSVERSE  SHIFT  VELOCITY  UPAR 
WRI TE (6*810) 

810  FORMAT(//,32X, ‘VALUES  FOR  UBAR‘*//> 

WRITE<6*820) 

ft  2C  FORMAT<3X*‘R**43X*‘  I  **/> 

WRI  TE (8,83  0  )  ( I  ST (I ) • 1  =  1  *  I  STD ) 

8  30  FOR  PAT (  4X  *6  I  12 ) 

no  850  j  =  t*  m 

WRI TE <6*840  >  R  <  J>  *  <UB(I*J)*I=1*ISTD> 

840  F  0  R  M  A  T ( F  7 • 3 *  6  E  1  2  •  4 ) 

850  CONTINUE 

URITE<6*852> 

852  FORMAT  <//*3X*‘  I  ‘*32X**R‘> 

WRI TE  <6  *  854 )  R ( 1 ) « R < I  Ml ) « R < I M 2 ) « R ( I  M3 ) • R C  I M4  ) 

854  F0RPAT<6X,SF12.4> 

DO  R58  I  =  1*  ISTO 

WRI TE(6,«)  I *UB< I«1)VUB(I«IM1), U9( I« IW2) »UB( I  *  I M3> 
t  *UB<1*IM4> 
ft  FOR  P  AT ( I5*5X*5E12»4) 

858  CONTINUE 

RESULTS  FOR  SECONDARY  FLOW  VELOCITY  U 
WRI TC<6*860 ) 

860  FORRAT(//»3?X*‘VALUES  FCR  U‘*//> 

URITE(6»8?0) 

WRITF(6«830)  (ISTU  >*I  =  1*ISTD> 

DO  ftftO  J  =  1  *  *1 

WRITE<6,840>  R(J>,  < U < I  * J > *  I  =  1  *  I  STD > 
per  CONTINUE 

RESULTS  FOR  STREAMWISE  VELOCITY  V 
URITE(6*890> 

8  S  0  F0RPAT<//,32X*»VALUES  FOR  V‘,//> 

WRI T E  C 6 *  8 2 3 > 

WRI TE(6*83C)  < I  ST < I > *  I =1  *  IS  TO > 

DO  9  0  0  J  =  l *  M 

WRITE<6,840>  P.  (  J  )  *  <  V  (  I «  J)  « 1  =  1 «  I  STD  > 

SOU  CONTINUE 

W»I TF(6» 852) 

WRITE  (6*854  )  P(l>*  R  (  I M 1  )  *  ft(IN?>»  P  <  I  *4  3  >  *  R(I*4) 

DO  ^CP  I  =  1*  ISTD 

WRITE(6»H)  I«V(I*1)*V(I*IM1)*V(I*I8I2)*V(I*IM3)*V(I*IM4) 
908  CONTINUE 


•  V*  *,*  -  •  ii  -  o  h  ' 
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C  RESULTS  FOR  LOCAL  FLOW  DEPTH 

URITE(f,,P10) 

9  10  F0RMAT(//»32X, ’VALUES  FOR  HL’,//> 

URI TE(6,820  ) 

WRITE(6,830)  ( I 5T( I ) »  1  =  1  ,  TSTQ) 

DO  920  J  =  1,  M 

WRITE (6,840)  R(J>,  < HL ( I * J ) « I  =  1 ♦ I  STD > 

<920  CONTINUE 

URITE(6»921 ) 

9  21  FORMAT!//, 25X, ’VALUES  FOR  UPAR  ♦  U*,//) 

WRI TE (6,820 ) 

URITE(6,3S0)  ( I  ST ( I > , I  =  1 , IS TD > 

00  930  J  =  l,  H 

WRITE (6,840)  R(J),  (UPPU( I, J) ,  1  =  1,  ISTP) 

9  3  0  CONTINUE 

URI TE(6,P31  ) 

931  F0°MAT(//,25X, ’VALUES  FOR  SQR T ( ( UB AR ♦U) * *?* V* *2 ) • « ft > 
WRITE (6, 820) 

WRITE(6»8304>  ( I  ST  ( I  )  ,  I  =  1 , 1 S  TD  ) 

DO  940  J=l,  M 

URITE(6,840)  R(J),  (UT*AN(I,J),  I=1,ISTD> 

9  4  G  CONTINUE 

WRITE(6,941 ) 

941  FORMAT!//, 25X, ’VALUES  FOR  VELOCITY  VECTOR  ANGLES’,//) 

WR I TE (fa , 820 ) 

URI TE (6, 83C )  C I  ST ( I ) , 1  =  1 , IS TD ) 

DO  'i 45  Jsl,  M 

WRI TE(6,840)  R(d),  ( ANGLE ( I , J ) ,  I=l,ISTD) 

945  CONTINUE 

URI TE(fc,946) 

946  FOR  HAT !//,2 OX, ’VALUES  FOR  UNIT  SEDIMENT  DISCHARGES’,//) 
URI TE(6,82n ) 

WRITE(6,83T)  ( I  S T ( 1  )  , I  =  1  ,  I  STD  ) 

DO  947  Jsl,  M 

URITE(6,84n)  R(J),  (QQSS(I,J),  I=1*ISTD) 

947  CONTINUE 
953  CONTINUE 

C 

C  PRIMOS  I/O  COMMENT 

C 

CLOSE(b) 

CLOSE (6) 

STOP 

END 

C  ••*»*****.•****»****•**»**,**.. 

SUBROUTINE  PQN(M,R,V,HL,VA,VQ,AT,QT) 

C  THIS  SU9R0UTTNE  DETERMINES  THE  FLOW  DISCHARGE  F^R  A  * 

C  GIVEN  RECTANGULAR  CROSS  SECTION  OF  CONSTANT  * 

C  BOTTOM  SLOP!.  WITH  KNOWN  VELOCITIES  AND  o P  T  H r  * 

C  VARIAHLES:  * 


K*. 
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r  **  r  NUMBER  OF  INCREMENTS  OR  VERTICALS  IM  THE  CKOSS  SECTION 
C  F  r  LOCAL  RADIAL  POSITION  ARRAY  (FT) 

r  V  =  LOCAL  FLOW  VELOCITY  ARRAY  (FT/SEC) 

C  Ft L  =  LOCAL  FLOW  DEPTH  ARRAY  (FT) 

C  VA  r  UNIT  FLOW  AREA  ARRAY  (St)  FT) 

C  VC  =  UNIT  FLOW  DISCHARGE  (CFS) 

C  AT  =  TOTAL  FLOW  AREA  (CU^IC  FT) 

C  QT  =  TOTAL  FLOW  DISCHARGE  (CFS) 

C****************************** 

DIMENSION  R (  1 7 ) •  V(17>*  HL(17)*  VA(17)*  VQ(17) 

AT  r  0.0 
QT  =  0.0 
DO  10  I  =  1 *M 
v  a ( r )  =  o.o 

VO (  I  )  r  0.0 
10  CONTINUE 

DO  ?  0  3  J  =  1  * M 
I F ( J  .EQ.  M)  GO  TO  100 
IF(HL(J+1)  .LT.  O.C)  GO  TO  120 
IF ( J  .GT.  1)  GO  TO  50 
W2  ;  0.5  *  ABS(  R(J*1)  -  R  ( J)  ) 

H2  =  (  HL ( J )  ♦  HL(J+1)  )  /  2.G 

V  A ( J )  =  w2  *  (  HL ( J )  ♦  H2  )  /  2.0 
VO<J)  =  VA(.J)  ■  V(J) 

GO  TO  no 

53  CONTINUE 
W 1  =  UP 
HI  =  H? 

W2  =  0.5  *  APS (  R ( J* 1 )  -  R(J)  ) 

H2  =  (  H L ( J )  ♦  HL(J-M)  )  /  2.0 

A1  =  U1  *  (  HI  ♦  HL ( J)  )  /  2.0 

A2  =  U2  *  (  HL(v')  *  H2  )  /  2.0 

V  A  (  J )  =  A  1  ♦  A2 
VQ  (  J )  =  V  A ( J )  *  V(J) 

CO  TO  110 

ICC  CONTINUE 
W1  -  hi  2 
hi  =  nr 

V  A ( J )  =  VI  *  (  HI  ♦  HL(J)  )  /  2.0 

V  Q ( J  >  =  V  A ( J  )  *  V(J) 

110  CONTINUE 

AT  =  AT  ♦  V  A ( J ) 

GT  =  QT  ♦  V Q  ( J ) 

GO  TC  2  n  0 
120  CONTINUE 
W1  =  W2 
HI  =  H? 

U 2  =  0.5  *  ABS  (  R ( J ♦  1  )  -  R ( J )  ) 

H2  -  (  HL ( J )  ♦  HL(JM)  )  /  2.0 
FAC  =  HI  t  (HI  -  HL  ( JO  )  ) 

WO  =  FAC  *  ACS(  R(JM)  -  (  R  ( J )  ♦  MJ-l)  )  /  2.C  ) 


I F  <  H2  •  LT .  0.0  )  GO  TO  150 
OELV  =  AOS (  UO  -  <  W1  ♦  U2  >  > 

VAC J*1 )  =  0.5  *  H2  *  DELW 
VOCJ*l>  =  VAC  J+ 1 )  *  VCJ) 

A1  :  U1  »  (  HI  ♦  HLCJ)  )  /  2.0 

A2  =  W2  *  C  HL  <  J )  ♦  H2  I  /  2.0 

VACJ)  =  A1  ♦  A? 

VO<J)  =  V A ( J )  *  V(J) 

AT  =  AT  ♦  V  A ( J )  ♦  VAC J*l) 

QT  =  QT  ♦  V3CJ)  ♦  V  G  C  J  ♦  1  ) 

GO  TC  2 20 
150  CONTINUE 

VACJ)  =  0.5  *  UO  *  HI 
VQC J)  =  VACJ)  *  VCJ) 

AT  =  AT  ♦  VACJ) 

OT  =  OT  ♦  VQCJ) 

GO  TO  220 
200  CONTINUE 
220  CONTINUE 
RETURN 
END 

C***«  ....  ******.*.*..***..*.  .*  *. 

SUB  ROUTINE  P6CFfVDARtD5  0*THETAC.POK.A*B*RHOS*G*Gl  *  r,2*G.\) 

C  THIS  SUBROUTINE  DETERMINES  DIMENSIONLESS  PARAMETERS  * 

C  Gl»  G 2 »  AMD  G3  FOR  THE  DEFINED  AND  INPUT  QUANTITIES  * 

C*.  **************** 

c  defined:  * 

C  RFCS  =  SPECIFIC  GRAVITY  OF  THE  SEDIMENT  PARTICLES  * 

C  B  =  PARAMETER  IN  THE  VELOCITY-SHEAR  RATIO  RELATION  * 

C  A  =  PARAMETFR  IN  THE  SHEAR-BED  LAYER  RATIO  RELATION  * 

C  INPUT  variables:  * 

C  F  =  D ARCY-WE ISB ACH  FRICTION  FACTOR  * 

C  VBAR  =  MEAN  STREAMWISE  FLOW  VELOCITY  CFT/SEC)  * 

C  D  50  =  MEDIAN  BED  MATERIAL  PARTICLE  SIZE  (MM)  * 

C  THE  TAC=  SHIELDS  CRITICAL  SHEAR  PARAMETER  * 

C  PCR  =  BED  LAYER  °GR0S1TY  * 

C  G  =  GRAVITATIONAL  CONSTANT  * 

C*************...***.******.'*.*. 
CALL  PG1 CF*B*G1) 

CALL  PG2CF,G2> 

CALL  PG3  CVBAR *RHOS*D50*F  »B t A  *  THE T AC  * POR t G * G3 > 

RETURN 

END 

C  *****  *  ******************  ******* 

SUBROUTINE  PG1CF  *BtGl  ) 

RN  =  1.0  /  SQRTCF) 

T  =  C  3  •  0  *RfJ  ♦  1.0)*C2.0*RN  ♦  1  •  0  )  /  C  2  •  0  *  R  N*  *  2  ♦  KM  ♦  1  .  P  ) 

Gl=T*B»F/fl.O 

RETURN 

END 


SIJ1'  R  CUT  I  NE  PG2(EtG^>) 

KN  =  1.0  /  S(JRT(F) 

T 1  r  (3.0  •  RN  ♦  l.C)  *  <2.0  *  R  N  ♦  1.0)  *  (  R  N  ♦  1.0) 

T  2  :  ( ?  •  0  •  RN *  * ?  ♦  R  N  ♦  1.0)  *  (RN  *2.0)  *  RN 
G?  r  T 1  /  T? 

RETURN 

END 

•  *  *  *************************** 

SUP* CUT  I NE  PG3  (  VBAR  »  R  HO  S  *D50*F*?*A«THET AC«P0R*G*G3  ) 

T1  =  F  *  THETAC 

T  2  =  S.O  *  G  *  050  *  (  R  HOS  -  1.3) 

T3  =  SGPT(T1  /  T?) 

T  A  =  r  *  VMAR  /  (A  *  C1.0  -  POR)) 

G  3  =  T  4  *  T  3 
return 

END 

******************  ************ 

SUBROUTINE  PF(F,H*RCtVBAP*Gl*C2*G3*R*ST*U*Fi#F2,F3*F4) 


THIS  SUBROUTINE  DETERMINES  ALL  FOUR  OF  THE  * 

GIMTN SIGNLESS  FUNCTIONS  FI*  F2*  F3*  ANO  FA  FOR  THE  * 
INPUT  QUANTITIES  GIVEN  * 

•  *  *  *************************** 
INPUT  variables:  * 

F  =  HARCY-UE ISBACH  FRICTION  FACTOR  * 

H  =  MEAN  DEPTH  OF  FLOW  OVER  CROSS  SFCTICN  (FT)  * 

RC  -  RADIUS  OF  CURVATURE  OF  CHANNEL  CENTERLINE  (FT) 

VBAR  =  MEAN  STR  E  AMU  I SE  FLOW  VELOCITY  (FT/SEC)  * 

G1 fO?»G3  =  DIMENSIONLESS  PARAMETERS  * 

R  =  RADIAL  POSITION  FROM  CHANNEL  CENTERLINE  (FT)  * 

ST  =  TRANSVERSE  BED  SLOPE  * 

J  =  MAXIMUM  SECONDARY  FLOJ  VELOCITY  (FT/SEC)  * 

F1.F2*R3*FA  =  DIMENSIONLESS  FUNCTIONS  TO  BE  EVALUATED  * 


***********  A******.*******.*** 

CALL  P  F 1 (ST  »RC*R  *H*F  1) 

CALL  PR? (ST*RC* R *GJ *G2*G3*H*F2) 

CALL  DP3(ST*RC,UfR»Gl*G2»G3,HfF,VBARtF3) 

CALL  PK A (ST «R«H« FA  ) 

R FT  URN 
E  NO 

**********  *  *******.*******  *  *4* 

SUBROUTINE  PF1 (ST *PC*R.H*F1 ) 

FI  =  ST  ♦  <  H  ♦  R  *  ST  )  /  (  RC  ♦  R  ) 

RETURN 

END 

***  ****  ************  ******  *  ***# 

SUBROUTINE  pF2(ST*RC*R*Gl*G.?*G3*H*F2) 

T 1  =  G2  *  G3  *  *  /  RC 

T?=R*G1*ST/H 

F2  =  RC  *  <  T1  -  T?  )  /  (  FT  ♦  R  ) 

R  ET URN 


r>  n  ci  n  r>  r>  o 


C********************  *******  *  *  *  * 

SUP ’OUT  INC  PF3(ST,RC,U,ft *G1 »G2*G3*H«F*VBAP,F3> 

T 1  =  F  *  U  /  (8.0  *  V  B  A  R  > 

T2  =  (  G2  *  H  )  /  <  RC  ♦  R  > 

T22  =  f,i  .  rc  *  ST  /  (  G3  *  (  RC  ♦  R  )  ) 

RN  =  1.0  /  SORT ( F ) 

T3=(R*ST*H)/H 
T 33  =  S3  *  R  *  U  /  (  H  *  VPAR  > 

T4  r  (  T2  -  T22  )  *  <  T2  -  T33  )  /  (  2.0  *  ^  ♦  1.0  > 

F  3  =  T  4  -  T  1 

RETURN 

END 

£*♦*****•**■**•*****•***********  **** 

SU3R  OUT  INC  PF4 ( ST  »R  ,H  *F4  ) 

T  =  R  *  ST  /  H 
F4  =  1.0  ♦  T 
RETURN 
E  NO 

c*************************.***** 

SUBROUTINE  EVAL (A*P  *C*RJ*RJM1 »  S  UM ) 

C  THIS  SUBROUTINE  EVALUATES  THE  ANALYTICAL  SOLUTION  OF  * 

THE  C  CONTINUITY  EQUATION  AS  THE  PROGRAM  BEGINS  A  NEW  * 

SECTION.  THE  GENERAL  FORM  OF  THE  INTEGRAL  IS  * 

C-IVEN  by:  * 

I  =  INTEGRAL!  X  *SGRT(  C*(A*X  ♦  P)/(X  ♦  C)  )  )  * 

where:  * 

X  =  RADIAL  COORDINATE *  R<J>  * 

A  =  TRANSVERSE  BED  SLOPE  *  ST  * 

3  =  MEAN  FLOW  DEPTH,  H  * 

C  C  =  RADIUS  OF  CURVATURE*  RC  * 

C  T  =  TRANSFORMED  COORDINATE*  DEFINED  AS:  * 

C  =  SORT!  C  *  (A  *  X  ♦  P)  /  (X  ♦  C)  >  * 

C  =  SORT!  RC  *  (  R ( J >  *  ST  ♦  H  >  /  (  R(U>  ♦  PC  )  >  * 

C  R J  *R  JM  l  =  UPPER  ANO  LOWER  INTEGRATION  LIMITS*  RESPECTIVELY  * 

C  SUM  =  FINAL  EVALUATION  OF  THE  INTEGRAL  * 

C***************************»*** 
T  2  =  SORT!  C*!A*RU*P)/(RJ  +  C)) 

T 1  =  SORT!  C  *  (  A  *  R JM 1  ♦  P  >  /  (  RJM1  ♦  C  )  > 

SUM  =0.0 
TT  =  A  *  C 

I F ! T  T  .NE.  0.0)  GO  TO  70 
WRITE (8*60)  A 

60  FORMAT!/, 5X* »ZrRO  TRANSVERSE  BED  SLOPE  ST  =  •*r20.x*/> 
STOP 

70  CONTINUE 

I F  <  T  T  .GT.  0.0)  GO  TO  100 

RIA  =  AT  AN!  T  2  /  SQRT(-TT)  )  -  AT  AN (  T1  /  SQPT(-TT)  ) 

RIA  =  RIA  /  SORT!  -TT  ) 

GO  TC  ISC 
100  CONTINUE 

RIA?  =  APS!  (  T  2  -  SORT!  TT  )  )  /  (  T2  ♦  SORT!  TT  )  )  ) 


B24 


RIA1  =  A8S(  (  T 1  -  SORT  (  TT>)/(  T 1  ♦  SrtPTt  TT  )  )  ) 
RIA  =  (  ALOSCRIA?)  -  ALOG(RIAl)  )  /  (  2.0  .  SOrJI 
l^O  CO*-l  T  I  MU  F  M,KT<  TT  )  > 

FAC  2  =  (  T?**2  -  TT  ) 

F  AC  1  =  (  T1  *  *2  -  IT  ) 

$1  :  1  12  /  FAC  2  **2  )  -  <  T l  /  F  AC1 *  *2  ) 

51  =  SI  *  (  tT  -  B  )  /  4.0 

52  =  (  T2  /  FAC2  )  -  {  T1  /  FAC1  ) 

S2  ~  S2  *  <  5.3  *  TT  -  B  )  /  {  8.0  *  TT  > 

*  J  ;*°  *  TT  +  p  >  *  RIA  /  (  8.0  *  TT  ) 

SUM  .  oO  *  J  TT  -  p  )  *  (  si  ♦  S2  ♦  S3  )  *  C**2 

R  LI  JRM 
f'NC 


OUTPUT  FILE:  OUTT 


MATHEMATICAL  MODEL  FQR  THE  PREDICTION  OF 
1HE  VELOCITY  FIELD  IN  RIVER  CLOW 


VALUES  FOR  VBAR,  H,  W,  SCL,  FOR,  3G,  RMU,  NSEG  : 

1.560  3. SOS  0  030  0  i 0400E-02  3.400  2.650  0.U0E-04  4 

FICTION  NUMBERS  WHERE  NEW  SEGMENTS  BEGIN  ARE  : 

1  137  273  40?  545 

uttttnmtttutnntttntuutttntntttnnHttttnttntntmnmttmnt 

DOWNSTREAM  STEPS  =  545  RADIAL  STEPS  =  17  CENTER  AT  M  =  ? 

RADIAL  POSITIONS  STORED  AT  J  =  1  5  9  13  17 

RESULTS  OUTPUT  EVERY  2  SECTIONS 

D/S  AND  RADIAL  OUTPUT  FREQUENCY  IS  1  1  STEPS 

RELATIVE  ERROR  CRITERIA  FOR  V  AND  UBAR  ARE  0.00100  0.01000 

MAXIMUM  ITERATIONS  20  PRINTED  EACH  10  ITERATIONS 

PROGRAM  OPTIONS  FOR  UBAR,  MOMENTUM  FORM,  DIRECTION  ARE  3  2  2 

INITIAL  NUMBER  OF  SUBINTESVALS  FOR  SIMPSON  RULE  =  0 

MAX  NUMBER  OF  U-V  ITERATIONS  IS  =  20 

DIMENSION  OR  NUMBER  OF  SECTIONS  TO  BE  TABULATED  IS  =  6 

AND  ARE  AT  SECTIONS. 

1  69  255  341  477  545 

SEDIMENT  POWER  LAW  OF  THE  FORM  Q3  =  A  *  (  V  )W 
WITH  A,  P  =  0.1080  40005 

DISCHARGE  =  6.3024  DARCY -WEISBACH  F  =  2.355483E-01 

POWER-LAW  N  =  0.4242?19flEl01  RADIAL  STEP  =  0  530DE+00 

LEFT  6  RIGHT  BANK  AT  R  =  0.4000E+01  -0.400CE+01 

SHEAR  VELOCITY  O.129971E+0C  FROUDE  NO  =  0 . 337C14E+0G 

N-TERM  GIVEN  BY  0.10377EE-01 
QS  =  D.639623E+00  ENG  TONS/DAY 

ttitnnmttiiititntntmmtttmnttmtnmmnmmnmntttmtm 


NHNDIMENSIONAL  QUANTITIES  . 

TWO  GRAVITY  TERMS  FOR  -VERT-,  HCR -  ARE  0  667648E+01  0  105766E+03 
DEPTH  =  0.631250E-01  LEFT  AND  RIGHT  BANKS  AT  0  500500E+00  -0  500000E+00 

RADIAL  STEP  0  625000E-0S 

ttttttttttittitttttttunttttttntttunnntnitttntttttntinttttnnntttt 


Ni  W  SEGMENT  —  IMPORTANT  PARAMETERS  GIVEN  AS 
FOR  SECTION  I  =  1  RC  =  3 . 433000E+02 


SEGMENT  LOCATION  BETWEEN  0  300000E+00  O.67SO0OE+O2 
ALPHA,  BETA  =  1.4160  3 . 2760  THETAC,  BSD  =  0  0320  0  3000 

NUMBER  OF  SECTIONS  BETWEEN  Si  4  52  IS  136 


COMPUTED  VARIABLES  FOR  THE  NEW  SEGMENT  GIVEN  BY 

RESTAR,  VST ARC,  RATIO  =  0  116313E+02  0.4989950-Oi  0.317931E+01 

DENSIMETRIC  FROUDE  =  0.632463E+01 

Gi ,  G2,  G3  =  O.7181779E-01  0  6249260E+00  0.3922577E+09 

tittitttntttttttttnmtttmntttntmtttmtxxutmmtttnmmimtn 

NONDIMENSIONALIZED  QUANTITIES  GIVEN  BY 


RC,  Si,  S2  =  3.5375050+01  O.OOOOOOE+OO  0.0437SOE+01  D50  0  194901E-02 

INTERVAL  BETWEEN  SECTIONS  DS  =  0  6204044E-01 

tniixxixxxtxxtttxtxxiuttxtxxxxxxxuxtxxtxumtxxtnixntxxixinzxxxxxxxxxtx 

INLET  SECTION  V  - 


0  10S0G30E+31  0.1943364C+0i  0.1036822E+01  0 . 103C402E+01  0 . 1024100E*01 
0.10179120+31  0.101183SE+01  3.i00586SE+Oi  0 .13003800+01  0.9942364E+00 
093857120+ 00  0.9830016E+00  5.9775252E+0O  0  9721394E+00  O.9663416E+P0 
0.96162950+00  3.95650030+03 


Q,  QT,  ETA  =  O.iOOOCOOOE+Ol  3.13013936E+91  0  99390757E+00 


ETA-MODIFIED  V-VELOCITIES  WIT!-;  UNIT  DISCHARGES 
J  =  1  V  =  0.10488832E+01  VA  =  0. 312500300-01  VQ  =  0.32777600E-01 

J=  2  V  =  0.10422237E+01  VA  =  0.625000000  01  VQ  =  0.6S138981E-O1 

J=  3  V  =  0  10356894E+01  VA  =  0. 625000000-01  VQ  =  0  6473C585E-01 

J  =  4  V  =  0.iC292764E+0i  VA  =  9.62500000E-01  VQ  =  0.64329773E-01 

J  =  SV=  9  10229809E+01  VA  =  0.625000000-01  VQ  =  0.63936308E-O1 
J  =  6  V  =  0.10167997E+01  VA  =  0 . 625008000-01  VQ  =  0.63549981E-01 

J  =  7  V  =  0.10107291E+01  VA  =  0 .625000000-01  VQ  =  0.631705670-01 

J  =  8  V  =  0.10047662E+01  VA  =  0 . 62500000E-01  VQ  =  0  62797389E-01 

J  =  9  V  =  0.99390757E+00  VA  =  0.625000000-01  VQ  =  0  62431723E-01 

J  *  10  V  =  0.99315012E+00  VA  =  0 . 62500000E- 01  VQ  =  0  62071882E-01 

J  =  11  V  =  0.98749113E+00  VA  =  O.625O0OOOE-O1  VQ  =  0  61718196E-01 

J  =  12  V  =  0  93192763E+D9  VA  =  0. 625000000-01  VQ  =  0 .61370477E-01 

J  =  13  V  =  0.77645724E+00  VA  =  0.62500000001  VQ  =  0  61023577E-01 

J  =  14  V  =  0.97107732E+00  VA  =  0. 625000000-01  VQ  =  0  60692333E-01 


J  =  15  V  =  C  9657S52<£+00  VA  = 

J  =  16  V  =  fl  960S7O92E+C0  VA  = 

J  =  1?  V  =  0  95545578E*flO  VA  = 

NEW  VALUE  Of  QT  WITH  MODIFIED  V  = 

95  =  0  64C432E>-00  ENG  TONS/DAY 


0  62508880E-81  VO  =  .0 .60361579E-01 
0.62500000C-01  VQ  =  0  60036182E-01 
0  31250000E  01  VQ  =  0  2?857??3F-0i 
0 . 1 OOOOOOOE  <01 

QS/QSACT  =  8  100440E+0I 


I  =  i  VALUES  FOP  UBAR  = 

0  000808DE+8C  0  3613371E-6 i 
0  12403 33E  ’■00  C  1343290E+00 
8  1283435E+03  0 . ii 68347C+00 
0 . 3530D59E  01  8  7L71790E-82 


0  66411S3E-01 
0 . 1392268E+0C 
8  101&711E+00 


0  91 11969E-01  0 . 11 050S1E+IO 
0  1397235E+CB  8. 1360422E+0C 
8  030241 IE-81  0  6104374E01 


END  OF  INLET  SECTION 


TOR  SECTION  I  =  3  ST  =  8  527754E-02  UCL  =  0  134S43E-01 

SATISFACTORY  ITERATION  ERRV2  =  0  226433E-03  ERRU2  =  8  417332E-02 

NUMBER  OF  ITERATIONS  KOUNT  =  2 

RELATIVE  DIFFERENCE  OF  9T  TO  S  IS  M61i8i5E-83 


CENTERLINE  POSITION  S  =  0  l240SlE>00  WITH  RC  =  8.S37S88E+01 

NEWEST  C  =  3  99998S46E+0B  WITH  ETA  =  0  1000461E+01 

SEDIMENT  DISCHARGE  =  0.64231F.+00  SS/3SACT  =  0  10042E+01 


NUMBER  OF  ITERATIONS  FOR  W,  0,  UV  ARE 
NEW  VALUES  FOR  UBAR  GIVEN  gr 


O.OOOOEiOG  0 .2017E-01  C.3617E -91  8.4863E-01 
C . 6716E-C1  8  6771E  Oi  fl .6S76E-81  0  6147E-01 
0  3606E-01  8.2336E-31  D.iOOiE  01  -3.S398E-82 
CIA- MODIFIED  MOMENTUM  VALUES  FOR  V  GIVEN  BY 
fl.i023E*Gi  0  1C44EV01  Q.l03CE*0i  0. 1031E-01 
0  10132+01  0  1007E+01  0  10C1E>01  0.9946E+0C 
0  9773E+00  0 . 9716E+-00  0 .?66iE+00  8  9606E+00 


0  S731E-C1 
0  550CE-O1 
-0  2222E-81 

0 . 1 025E+01 
0 . 9808E-*  03 
0  9552E+00 


0  6392E-01 
0  4648E-01 


0. 1019E+01 
0  9330E+80 


FOR  SECTION  I  =  S  ST  =  0.98B028E-02  UCL  =  8  2S1373E-01 

SATISFACTORY  ITERATION  ERRV2  =  0 . 191607E-03  ERRU2  =  0.489446E-02 

NUMBER  OF  ITERATIONS  KOUNT  =  2 

RELATIVE  DIFFERENCE  OF  9T  TO  G  IS  8  4338026E-03 

CENTERLINE  POSITION  S  =  0  24C162E+00  WITH  RC  =  0 . S37S00E+O1 

NEWEST  8  =  0 . 99998736E+00  WITH  ETA  =  0  1000434E+01 

SEDIMENT  DISCHARGE  =  8  64227E+00  QS/QSACT  =  0  10041E+01 


NUMBER  OF  ITERATIONS  FOR  UB.  Q,  UV  ARE 


o 


3 


IO  VALUES  FOR  UBAR  GIVEN  BY 

8.8088E+80  S . 13Q2E-C1  8.3331E-81  8.4443E-81  9.S249E-01  8.S772E-81 

8.683SE-01  8.68S7E-81  3.S8S7E-B1  I.54S3E-01  0.4862E-81  8  4897E-01 

8  3174E-81  0  2186E-91  -8.9C33E-82  -8.4213E-82  -O.iBSBE-Oi 
ETA-NODIFIED  HONENTUN  VALUES  FOR  V  GIVEN  BY 
8. i8B8£+0i  8.1044E+81  8.i039E+Bi  6.iC33E+0i  0  1027E<01  0  iC20E+Ci 

B.iOiAEm  8  1C93E+81  0.i302E*8i  8  996C£*B3  8.9933EM)8  8.984iE«-00 

8  9733E+00  8.9725E+80  G.9663E+80  0  96iiE»00  8  9S55E+83 


r  / 


VALUES  FOP  UBAR 


R  ) 

i  6?  200  341  47?  54S 

0  GOC  0  00900*30  800000+03  3  OOCCE+OO  0  1182E-0Z  0  1737E-03  8  6O42E-04 

0430  0  36140-31  0.1441E-02  8.2A29E-Q3  -fl  82400-03  -0  1618E-03  -O.6320E-I4 

0  370  0  66410-01  0  2746E-02  3.4+650  03  -8.26210-02  -0  4552E-03  -0  1663E-03 

8313  8  ?! 12E  11  0  37080-02  8  60600-93  -8  42090-32  -0  7049E-03  -9. 2498E-03 

0  200  0  HOSE i 3?  0.42590-02  3. 70920-03  -0  0S990-32  -O.913SE-03  -0  3150E-03 

0  100  0  1243080  3  47S4E  -02  0  St??E-03  -0  6771E  0?  -0  10820-82  -0  3623E-03 

0.120  8  13430*35  0.49340-02  8  V7060  03- -9  77340  32  -8  12090-82  -8  3931E-03 

0  043  8  1392080  8  49230  02  3  96220-83  -0.80630  82  -  0  12942-02  -0  4364E-03 

3880  313970 ‘00  8.47600-02  S.9S21E-03  -8  91122-02  -0  13300-02  -O.4831E-03 

8.863  8  13600*00  O.4449E-02  0  909OE-O3  -9  94940-02  -0.1327E-G2  -0.383OE-I3 

3.120  8  12330*08  0.48120-32  O  93080-03  -0  9406E-02  -O  1271E-02  -C  3467E-03 

8.138  0  1160E- 30  3.3464E-02  0.73020-03  -3  90720-82  -0  11600-02  -C.2963E-G3 

8208  0  10170  83  8  23i7E-02  3. ofl9flr  03  -0.83350-82  -0  7956E-C3  -fl  2344E-J3 

3.313  3  33820  01  0  23320-02  3. 4094E-03  -0.71310  02  -0  7S22E  03  -0.1653E-03 

8  370  9.61040-81  0  12690-92  8  23Q0E-D3  -3.0207E  02  -0  0311E  83  -D.9062E-J4 
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APPENDIX  C:  NOTATION 


NOTATION 

Constant  used  in  sediment-transport  formula 

Constant  in  quadratic  equation 

Constant  used  in  sediment-transport  formula 

Constant  in  quadratic  equation 

Constant  in  quadritic  equation 

Local  flow  depth 

Centerline  flow  depth 

Median  bed-material  size 

Darcy -Wei sbach  friction  factor 

Arbitrary  variable 

Dimensionless  functions 

Froude  number 

Gravitational  constant 

Dimensionless  functions 

Local  bed  elevation 

Centerline  bed  elevation 
Local  water-surface  elevation 
Streamwise  grid  locations 
Streamwise  index  for  section  numbers 
Transverse  grid  locations 
Transverse  index  for  section  numbers 
Number  of  transverse  grid  points 

Exponential  parameter  used  in  power-law  velocity  distributi 


ts 
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Number  of  streamwise  grid  points 
Bed-material  porosity 
Total-load  discharge  per  unit  width 
Total -load  discharge 
Total  water  discharge 
Transverse  coordinate 
Radius  of  curvature  at  inside  bank 
Radius  of  curvature  at  outside  bank 
Radial  grid  locations 
Centerline  radius  of  curvature 
Streamwise  coordinate 
Origin  of  streamwise  coordinate 
Streamwise  grid  locations  along  channel  centerline 
Streamwise  water-surface  slope  along  channel  centerline 
Transverse  bed  slope 
Integral  function 
..  Integral  functions 

Local  secondary-flow  velocity 
Local  shear  velocity 

Local  secondary-flow  velocity  at  water  surface 

Secondary -flow  velocity  at  channel  centerline 

Mass-shift  velocity 

Local  streamwise  velocity 

Depth -averaged  streamwise  velocity 

Depth-averaged  velocity  at  channel  centerline 
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T  Area-averaged  streamwise  velocity 

W  Channel  width 

y  Vertical  coordinate 

yb  Bed-layer  thickness 

a  Proportionality  constant  for  bed-layer  thickness 

$  Proportionality  constant  for  bed-shear  stresses 

9  Critical  shear-stress  parameter 

p  Fluid  mass  density 

ps  Sediment  mass  density 

t  Streamwise  bed-shear  stress  =  t 

o  os 

t  Transverse  bed-shear  stress 

or 

4  Streamwise  angular  coordinate 
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